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Abstract 


Detonation-based  combustors  leverage  the  higher  thermodynamic  efficiency  of  the  Atkinson 
cycle  compared  to  the  traditional  deflagration-based  combustion  of  the  Brayton  cycle.  The 
rotating  detonation  engine  (RDE)  has  one  or  more  shock  waves  rotating  around  an  annulus. 
The  RDE  can  theoretically  be  20%  more  thermally  efficient  than  a  traditional  deflagration- 
based  cycle.  An  RDE  was  modeled  in  Numerical  Propulsion  System  Simulation  (NPSS) 
based  on  a  model  developed  in  Microsoft  Excel.  The  thermodynamic  analysis  of  the  RDE 
in  these  models  is  broken  into  four  streams.  Empirical  models  were  used  to  fold  the  per¬ 
centage  of  the  total  flow  in  each  stream.  The  pre-detonation  pressure  was  iterated  until 
the  entrance  mass  flow  calculations  matched  the  exit  mass  flow  calculations.  A  parametric 
analysis  was  used  to  compare  the  variation  in  specific  impulse  from  the  NPSS  model  to  the 
Microsoft  Excel  model  and  other  published  results.  The  RDE  has  a  peak  air-breathing  en¬ 
gine  specific  impulse  of  approximately  5,500  sec  and  a  peak  rocket  engine  specific  impulse 
of  approximately  150  sec. 
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COMPUTER  MODELING  OF  A  ROTATING  DETONATION  ENGINE  IN  A  ROCKET 


CONFIGURATION 

I.  Introduction 

Conventional  jet  and  rocket  engines  extract  energy  from  fuel  using  deflagration  processes, 
which  follow  a  Brayton  cycle  [6].  The  efficiency  of  an  ideal  Brayton  cycle  depends  on  the 
ratio  of  temperatures  at  the  entrance  and  exit  of  the  process  [6],  so  increasing  the  exit 
temperature  increases  the  efficiency  of  the  process.  Material  properties  limit  the  maximum 
exhaust  temperature,  thus  limiting  efficiency  [7]. 

A  detonation  cycle,  however,  follows  an  Atkinson  cycle  [8].  Figure  1  illustrates  the 
thermodynamic  nature  of  the  Brayton  and  Atkinson  cycles  via  the  temperature-entropy  (T- 
S)  and  pressure-specific  volume  (P-is)  diagrams.  The  superior  performance  of  the  Atkinson 
is  evidenced  by  the  greater  area  under  the  two  curves  in  Fig.  1.  For  equal  heat  addition  and 
mechanical  compression,  the  Atkinson  cycle  is  more  thermodynamically  efficient  than  the 
Brayton  cycle  [8]. 

One  utilization  of  a  detonation  cycle  is  the  rotating  detonation  engine  (RDE).  Research 
into  the  RDE  began  as  early  as  the  1950s  at  the  Lavrentyev  Institute  of  Hydrodynamics  in 
Siberia  [9].  A  major  motivator  for  the  study  of  RDEs  is  their  increased  thermal  efficiency; 
the  RDE  can  theoretically  be  20%  more  efficient  than  a  traditional  deflagration-based  cycle 
[2;  10;  11]. 

In  an  RDE,  one  or  more  shock  waves  rotate  around  an  annulus.  Figure  2  illustrates  a 
temperature  contour  plot  from  a  computational  fluid  dynamics  (CFD)  flowfield  visualization 
of  an  RDE  simulation.  Fuel  and  oxidizer  entering  from  the  inlet  is  detonated  at  the  shock 
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Figure  1.  Comparison  of  thermodynamic  diagrams  for  Brayton  and  Atkinson  cycles  [1] 


wave  traveling  around  the  annulus  and  is  expelled  out  of  the  exhaust  end. 

Since  the  radial  thickness  of  the  RDE  is  significantly  smaller  than  the  length  or  circum¬ 
ferential  dimensions,  radial  variation  is  negligible  [2-5;  10-14].  The  three-dimensional  RDE 
in  Fig.  2  can  be  “unrolled”  into  the  two-dimensional  RDE  in  Fig.  3. 


Exhaust 


Figure  2.  Temperature  contour  plot  of  an  RDE  simulation  [2:  Figure  1] 
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Figure  3.  Two-dimensional  “unwrapped”  temperature  contour  plot  of  the  RDE  [3:  Fig.  2] 
Problem 

Kaemming  developed  a  preliminary  model  of  an  RDE  using  Microsoft  Excel  and  Microsoft 
Visual  Basic  macros  [4],  The  objective  of  this  thesis  was  to  develop  a  similar  RDE  model  in 
Numerical  Propulsion  System  Simulation  (NPSS)  and  describe  the  utilization  of  the  NPSS 
model  in  a  rocket  configuration.  The  rocket  configuration  consists  of  a  flow  inlet,  the  RDE, 
and  the  exhaust. 

Methodology 

The  Erst  step  in  the  analysis  of  the  RDE  is  the  translation  of  the  RDE  element  based 
on  the  Kaemming  model  [4]  into  an  NPSS  element.  The  new  NPSS  element  is  incorporated 
into  a  full  RDE  model  by  initializing  and  attaching  the  air  inlet,  fuel  inlet,  and  exhaust  fluid 
ports.  The  results  of  running  the  NPSS  model  were  compared  to  results  from  the  Kaemming 
model  [4;  15]  and  published  CFD  simulations  [3;  14], 


3 


II.  Previous  Work  on  Rotating  Detonation  Engines 


This  section  describes  the  station  numbering,  a  thermodynamic  overview  of  the  RDE, 
and  the  empirical  models  used  in  the  RDE  model. 


Station  Designation 


A  station  designation  consistent  with  previously  accepted  station  numbering  designations 
has  been  agreed  upon  by  Joint  Army  Navy  National  Aeronautics  and  Space  Administra¬ 
tion  (NASA)  Air  Force  (JANNAF)  (in  JANNAF-GL-2013-01)  and  Society  of  Automotive 
Engineers  (SAE)  (in  SAE  AS755f)  [15;  16],  illustrated  in  Fig.  4.  In  the  rocket  configuration 
implemented  in  this  thesis,  the  inlet  exit,  air  valve  exit,  and  RDE  entrance  are  assumed  to 
be  collocated  (i.e.,  station  3  is  the  same  as  3*). 
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Figure  4.  Station  designations  used  in  rocket  configuration^:  Fig.  3-1] 
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Thermodynamic  Overview 


There  has  been  considerable  thermodynamic  analysis  of  the  RDE  [2;  5;  11;  13],  so  only  a 
summary  is  provided  here.  The  main  process  taking  place  in  the  RDE  is  the  detonation  wave. 
A  steady-state,  one-directional  planar  detonation  can  be  modeled  based  on  Zel’dovich,  von 
Neumann,  and  Doring  (ZND)  theory  [2;  5;  11;  13].  In  the  “unrolled”  RDE,  the  detonation 
wave  is  a  traveling  planar  wave,  so  the  traditional  ZND  model  does  not  apply  in  its  present 
state  [2;  11;  13]. 

Nordeen  [11]  modified  the  ZND  model  using  velocity  triangles  to  change  the  frame  of 
reference  from  the  laboratory  frame-of-reference  to  the  detonation  wave  frame-of-reference. 
A  graphical  representation  of  a  velocity  triangle  is  illustrated  in  Fig.  5.  In  Fig.  5,  Uwave  is 
the  azimuthal  wave  speed  traveling  to  the  right  in  the  laboratory  frame-of-reference  and  V  is 
the  time-averaged  velocity  of  the  main  flow  [5].  The  normal  velocities  of  the  flow  upstream 
and  downstream  of  the  detonation  front  are  given  by  W\  and  W2 ,  respectively  [5].  In  Fig.  5, 
D  is  the  normal  speed  of  the  detonation  wave  in  the  laboratory  frame-of-reference,  and  is 
equal  and  opposite  to  W\  [5].  In  the  detonation  wave  frame-of-reference,  the  flow  appears 
steady  and  the  ZND  theory  can  be  applied  [5;  11]. 


Figure  5.  Velocity  triangles  used  to  change  between  laboratory  frame-of-reference  and  deto¬ 
nation  wave  frame-of-reference  [5:  Fig.  5] 

Schwer  and  Kailasanath  tracked  the  paths  of  fluid  particles  through  a  CFD  simulation 
of  the  RDE  in  the  detonation  frame-of-reference  [3].  Schwer  and  Kailasanath  defined  three 
distinct  thermal  regions  in  the  solution:  Detonation-A  (detonation  followed  by  an  oblique 
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shock  wave),  Detonation-B  (detonation  only),  and  Non-detonation  (deflagration)  [3].  The 
pathlines  and  three  thermal  regions  are  illustrated  in  Fig.  6. 


Figure  6.  Two-dimensional  “unrolled”  RDE  with  tracked  pathlines  and  three  thermal  regions 
[3:  Fig.  8] 

Kaemming  extended  the  results  from  Schwer  and  Kailasanath  [3]  with  the  introduction 
of  a  mixed  region  [4],  Additionally,  Kaemming  added  substation  nomenclature  to  organize 
each  of  the  steps  along  each  of  the  regions.  Kaemming’s  [4]  nomenclature  is  illustrated  in 
Fig.  7.  The  letters  a-d  describe  the  paths  through  the  thermal  regions:  path  a  has  only 
detonation,  path  6  has  detonation  followed  by  an  oblique  shock,  path  c  has  deflagration,  and 
path  d  has  a  mixture  of  paths  a,  6,  and  c  [4], 

The  air  and  fuel  enter  the  RDE  at  station  3*  and  are  injected  into  the  flow  to  substation 
3.2.  Substation  3.2  properties  apply  to  all  flow  prior  to  combustion,  and  all  of  the  paths 
pass  through  it  [4],  Table  4  summarizes  the  processes  that  take  place  between  substations. 
The  flow  traveling  along  paths  a  and  b  pass  through  a  detonation  between  substation  3.2 
and  a3.4  and  63.4,  respectively  [4],  The  flow  traveling  along  path  a  continues  to  travel 
to  the  exit  of  the  RDE.  The  flow  traveling  along  path  b  passes  through  an  oblique  shock 
between  substations  63.4  and  63.6,  and  then  continues  to  travel  to  the  exit  of  the  RDE  [4], 
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The  flow  traveling  along  path  c  is  deflagrated  between  substations  3.2  and  c3.4,  and  then 
passes  through  an  oblique  shock  between  substations  c3.4  and  c3.6[4].  All  of  the  flow  at 
substation  c3.6  mixes  with  portions  of  the  flow  from  substations  a3.4  and  63.6;  the  mixed 
flow  is  substation  <73.6.  The  percentage  of  the  flows  from  substations  a3.4  and  63.6  entering 
the  mixed  region  is  dictated  by  empirical  models  [4].  The  three  remaining  paths  a3.4,  63.6, 
and  <73.6  travel  to  their  respective  substation  exits:  aA,  64,  and  <74.  The  RDE  exit  station  4 
is  a  combination  of  the  substation  exits  a4,  64,  and  dA. 

The  analysis  of  each  path  is  described  in  Chapter  III. 
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Table  4.  Description  of  particle  movement  between  substations 


From 

To 

Process 

3* 

3.2 

Air  and  fuel  injected  into  RDE 

3.2 

<2:3.4 

Detonation 

63.4 

Detonation 

c3.4 

Deflagration 

63.4 

63.6 

Oblique  shock 

c3.4 

c3.6 

Oblique  shock 

a3.4  (partial) 

63.6  (partial) 
c3.6  (all) 

d3.6 

Mixing  based  on  empirical  models 

a3.4  (remaining) 

a4 

Exit  RDE 

63.6  (remaining) 

64 

Exit  RDE 

d3.6  (all) 

dA 

Exit  RDE 

Empirical  Models 

Kaemming  developed  several  empirical  models  based  on  CFD  simulations  to  reduce  the 
overall  complexity  of  the  solution  procedure  [3;  4;  14].  The  models  are  discussed  below. 


Plenum  Model. 

CFD  simulations  [3;  14]  show  that  the  rotation  of  the  detonation  wave  around  the 

annulus  of  the  RDE  causes  a  variation  of  static  pressure  P3.2  at  the  air  and  fuel  injectors 
[3;  14].  The  thermodynamic  properties  of  the  remainder  of  the  RDE  depend  on  P32,  so  a 
model  was  required  to  account  for  this  variation  along  the  plenum.  Kaemming  developed 
the  empirical  model  seen  in  Eq.  1  and  Eq.  2  for  the  injector  pressure  profile  based  on  the 
results  from  CFD  simulations  [3;  4;  14].  The  value  of  P32  varies  with  azimuthal  location, 
but  is  based  on  the  maximum  value  of  P3.2(0),  which  is  the  value  of  P3  2  used  in  the  current 
iteration  (P3.2  is  further  described  in  Chapter  III).  In  Eq.  1,  the  term  6  is  the  fractional 
pressure  drop  and  Tdr0p  is  the  time  until  the  first  drop  [4],  In  Eq.  2,  the  term  (p^)  is  the 


ratio  of  static  pressures  after  and  before  the  detonation. 


k=-'n{1-b)  (!) 

T dr  op 

Ps.2(t)  =  iVa(O)  |l+  (^)  -1  e-“|  (2) 

Published  results  [3]  showed  the  pressure  prohle  as  a  function  of  azimuthal  location,  while 
the  empirical  model  in  Eq.  1  and  Eq.  2  showed  pressure  prohle  as  a  function  of  time  [4], 
Kaemming  used  the  detonation  velocity  to  convert  the  time  to  distance,  and  the  RDE  cir¬ 
cumference  to  convert  the  azimuthal  distance  to  azimuthal  location  [15].  The  final  prohle 
was  translated  by  180°  and  the  pressure  was  converted  from  pounds  per  square  inch  to  atmo¬ 
spheres  to  match  the  CFD  results.  Kaemming  found  that  b  =  0.80  and  Tdrop  =  0.00005  sec 
produced  a  pressure  prohle  most  closely  matching  CFD  results  [4], 

The  empirical  models  from  Eq.  1  and  Eq.  2  are  used  in  the  NPSS  model.  A  comparison 
of  the  pressure  prohles  at  station  3.2  in  the  Kaemming  and  NPSS  models  to  the  prohle  found 
from  the  CFD  simulations  is  illustrated  in  Fig.  8.  The  pressure  prohle  example  in  Fig.  8  is 
based  on  a  specihc  case  using  the  same  input  parameters  as  CFD  simulations,  summarized  in 
Table  6.  Overall,  the  simple  model  follows  the  CFD  results  closely;  however,  two  regions  of 
significant  deviation  can  be  seen  at  azimuthal  locations  of  85°  and  165°.  Kaemming  claims 
the  deviations  are  the  result  of  second-order  effects,  such  as  secondary  waves,  that  may  affect 
operation  of  a  real  RDE,  but  do  not  affect  performance  analysis  [4],  The  agreement  of  the 
Kaemming  model  with  CFD  simulations  support  the  claim  [3;  14], 


Detonation  Height  Ratio  Empirical  Models. 

Kaemming  developed  several  empirical  models  that  are  functions  of  the  ratio  of  the  cal¬ 
culated  detonation  height  to  the  user-defined  RDE  height  ^  h^E  )  Kaemming  developed 
the  empirical  models  by  simplifying  CFD  solutions  and  prohles  [4;  15].  The  non-dimensional 
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Figure  8.  Comparison  of  pressure  profile  from  CFD  results  to  model  from  Eq.  1  and  Eq.  2 
[4:  Fig.  6-1;  3:  Fig.  9] 


ratio  of  heights  allows  the  RDE  model  to  be  scaled;  however,  experimental  comparisons  have 
not  been  accomplished  to  determine  the  accuracy  of  the  scaling. 

In  the  NPSS  model,  the  empirical  models  are  implemented  in  the  form  of  tables,  where 
the  program  calculates  the  value  of  and  interpolates  within  the  models  shown  in  Fig.  9. 


The  “Detonation  Plus  Shock  Flow  Fraction”  model  in  Fig.  9  provides  the  fraction  of  the 
total  detonated  flow  entering  the  oblique  shock  [4],  The  model  is  mathematically  described 
in  Eq.  3. 

Detonation  Plus  Shock  Flow  Fraction  =  - -  pa,tl1.6 -  (3) 

dlpath  a  T  blpath  b 

The  value  returned  from  Fig.  9  provides  the  percentage  of  detonated  flow  that  follows  path  b. 
The  remainder  of  the  flow  follows  path  a. 

The  “Relative  Mixing”  model  in  Fig.  9  provides  the  relative  amount  of  detonated  flow 
that  mixes  with  the  deflagrated  flow  [4].  The  model  is  mathematically  described  in  Eq.  4. 

Relative  Mixing  _  ^detonated  flow  mixing  with  path  d  /^\ 

dlpath  c 
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Height  Detonation , 

7 Height  RDE 

Figure  9.  Empirical  models  based  on  the  ratio  of  detonation  height  to  RDE  height 


The  value  returned  from  Fig.  9  provides  the  ratio  of  mixed  detonated  flow  to  total  deflagrated 
flow  from  path  c.  Kaemming  approximates  that  half  of  the  detonated  flow  entering  the 
mixing  region  is  from  path  a  and  the  other  half  from  path  b  [4].  All  of  the  deflagrated  flow 
from  path  c  enters  the  mixing  region,  so  the  total  mass  flow  entering  path  d  is  given  by  Eq.  5. 

?hpath  d  =  bldetonated  flow  mixing  into  path  d  "F  blpath  c  =  TTlp&th  c  (1  T  Relative  Mixing)  (5) 


The  “Circumferential  Velocity  -  Detonated  Flow”  and  “Circumferential  Velocity  -  Mixed 
Flow”  models  in  Fig.  9  provide  the  mass-averaged  circumferential  (azimuthal)  velocity  for 
the  detonated  and  mixed  flows  relative  to  the  total  velocity  of  the  detonated  or  mixed  flows, 
respectively.  The  “Circumferential  Velocity  -  Detonated  Flow”  and  “Circumferential  Velocity 
-  Mixed  Flow”  models  are  mathematically  described  in  Eq.  6  and  Eq.  7,  respectively. 

r  j'lTri  •.  ~rx  .  ,  ,  -r-,,  detonated, circumferential  /n\ 

Circumferential  Velocity  -  Detonated  1  low  =  - - : - - -  (o) 

detonated, axial  T  detonated, circumferential 
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/-v  r  i  •  i  T7-  i  •.  tv /r*  l  T~n  mixed, cir  cum  f  erential  /^\ 

Circumferential  Velocity  -  Mixed  r  low  =  — - : - - -  (7) 

'W'mixed, axial  “1“  ^ mixed, circumferential 

The  percentage  of  circumferential  velocity  for  the  detonated  flow  increases  with  non-dimensionalized 
detonation  height,  as  seen  in  Fig.  9.  A  higher  proportion  of  h^E  means  the  detonation  takes 
up  more  of  the  RDE  chamber,  and,  conversely,  the  non-detonated  flow  takes  up  less  of  the 
RDE  chamber.  The  shock  waves  try  to  align  the  exiting  flows,  and  the  smaller  non-detonation 
region  causes  a  greater  percentage  of  the  flow  to  travel  in  the  circumferential  direction  [4]. 
Conversely,  with  the  mixed  flow,  the  percentage  of  rotating  flow  decreases  as  the  detonation 
height  increases  because  the  flow  mixes  in  a  smaller  region,  requiring  less  rotation  [4], 

The  flow  exiting  the  RDE  is  highly  irregular  [3;  4;  14;  17].  The  “Velocity  Distortion” 
model  in  Fig.  9  characterizes  the  exit  velocity  using  the  maximum,  minimum,  and  mass- 
averaged  velocities  mathematically  described  in  Eq.  8  [4],  The  Kaemming  and  NPSS  models 
both  use  a  constant  velocity  distortion  of  120%  [4]  based  on  exit  velocity  profiles  found  from 
CFD  simulations  [14], 

Velocity  Distortion  =  ^ max  — Yum  (3) 

*  avg 

Kaemming  also  found  that  the  exit  mass  flows  in  the  CFD  solutions  [3;  14]  have  the  same 
corrected  mass  flow  [15].  In  the  Kaemming  model,  the  exit  flow  is  assumed  to  be  choked 
[4],  However,  the  choked  flow  assumption  calculates  an  exit  mass  flow  and  momentum 
flux  higher  than  the  actual  values  because  of  the  velocity  distortion  [15].  Kaemming  uses 
correction  terms  for  the  mass  (^pr)  and  momentum  (pqpr)  based  on  variations  in  the  velocity 
distortion  [4],  The  mass  and  momentum  correction  terms  are  mathematically  described  in 
Eq.  9  and  Eq.  10,  respectively,  and  are  illustrated  in  Fig.  10  [4],  To  find  the  actual  exit  flow, 
the  correction  terms  are  applied  to  the  calculations  where  choked  exit  flow  was  assumed  [15] . 

For  the  constant  velocity  distortion  of  120%,  the  correction  terms  are  treated  as  constants: 
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m 

m* 


0.8682. 


=  0.8555  and 


—  =  [—0.1004  (Velocity  Distortion)  +  0.0001]  (Velocity  Distortion)  +  1  (9) 

m* 

TTlV 

— —  =  [—0.0918  (Velocity  Distortion)  +  0.0003]  (Velocity  Distortion)  +  1  (10) 

1 1  i  V 


Velocity  Distortion 


Figure  10.  Correction  factor  for  the  mass  (m/m*)  and  momentum  (mV/mV*) 


Summary 

A  station  designation  has  been  agreed  upon  by  JANNAF  and  SAE  that  conforms  to 
traditional  engine  station  numbering.  The  thermodynamic  analysis  of  all  of  the  streamlines 
in  the  RDE  is  reduced  to  four  paths.  Several  empirical  models  are  used  to  simplify  the 
complex  flow  interactions  in  the  RDE. 
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III.  Model  Summary 


The  NPSS  RDE  model  is  composed  of  four  elements:  the  air  source,  the  fuel  source, 
the  exhaust,  and  the  RDE  element  itself.  The  flow  into  the  RDE  includes  the  air  and  fuel 
elements,  while  the  exhaust  element  receives  all  of  the  exhaust  from  the  RDE  element.  The 
air  source,  fuel  source,  and  exhaust  elements  are  all  ideal  components  with  no  losses,  as  are 
all  the  ports  that  link  these  elements  to  the  RDE  element.  In  all  cases,  hydrogen  fuel  (H2) 
is  used,  though  the  current  model  allows  for  the  use  of  ethylene  ( C2H4 ). 

The  functions  utilized  to  perform  repeated  calculations  are  summarized  in  Table  5.  De¬ 
tailed  descriptions  of  these  functions  can  be  found  in  Appendix  A. 

Table  5.  Summary  of  functions  implemented  in  RDE  model 


Function 

Description 

PtqP 

Isentropic  pressure  ratio 

TtqT 

Isentropic  temperature  ratio 

M_PtqP 

Calculate  Mach  number  from  isentropic  pressure  ratio 

M.TtqT 

Calculate  Mach  number  from  isentropic  temperature  ratio 

AqAstar 

Isentropic  area  ratio  for  choked  flow  in  a  nozzle 

Msub_A 

Calculate  subsonic  Mach  number  from  isentropic  area  ratio 

Msup_A 

Calculate  supersonic  Mach  number  from  isentropic  area  ratio 

P2qPlns 

Calculate  pressure  ratio  across  normal  shock 

Pt2qPtlns 

Calculate  total  pressure  ratio  across  normal  shock 

T2qTlns 

Calculate  temperature  ratio  across  normal  shock 

M2ns 

Calculate  Mach  number  across  normal  shock 

f ind_f low 

Find  mass  flow  to  satisfy  total  area  requirement  for  three  separate  streams 

expans i on_mode 

Find  exit  Mach  number  for  expanding  flow  from  throat  to  exit 

CJ_mach_area 

Calculate  Chapman- Jouguet  Mach  number  with  area  relief 

OneDNtrp 

Interpolate  along  arrays 

RDE  Element  Summary 


The  Kaemming  model  is  divided  into  three  solution  stages:  entrance  flow,  thermodynamic 
analysis,  and  exit  flow  [4],  The  NPSS  RDE  element  closely  follows  the  methodology  of  the 
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Kaemming  model.  The  individual  total  pressure  of  the  air  (Ptz^air)  and  fuel  (Ptzjuei)  entering 
the  RDE  element  are  known,  but  the  pressure  of  the  air  and  fuel  mixture  injected  into  the 
RDE  (i.e.,  Pt 3)  is  unknown.  The  RDE  element  iterates  P3.2  until  the  calculated  mass  flow 
entering  the  RDE  matches  the  calculated  mass  flow  exiting  the  RDE.  Figure  11  summarizes 
the  implementation  of  the  Kaemming  model. 


Station  Station 


The  default  parameters  and  values  for  user-defined  inputs  to  the  RDE  element  are  sum¬ 
marized  in  Table  6.  Table  6  also  summarizes  the  mathematical  and  NPSS  variable  names  of 
the  parameters.  These  parameter  values  must  be  changed  in  the  element  code  itself,  not  the 
run  hie.  The  values  for  the  parameters  in  Table  6  are  based  on  RDE  geometry,  model  test 
conditions,  or  from  Kaemming’s  analysis  of  CFD  simulations  [3;  4;  14]. 

Since  the  how  direction  is  assumed  to  be  two-dimensional,  an  ^-coordinate  system  is 
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used.  The  x-direction  is  the  circumferential  or  azimuthal  direction,  while  the  ^/-direction  is 
the  axial  direction  from  the  inlet  to  exhaust.  In  several  cases,  different  frames-of- reference  are 
implemented.  Static  thermodynamic  properties  are  independent  of  the  frame-of-reference, 
while  total  thermodynamic  properties  and  relative  velocities  are  affected  by  the  frame-of- 
reference.  Variables  used  for  analysis  in  the  laboratory  frame-of-reference  do  not  have  any 
additional  subscripts.  Variables  used  for  analysis  in  the  detonation  wave  frame-of-reference 
have  the  subscript  w. 
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Table  6.  Default  values  for  user-defined  input  parameters 


Parameter  Description 

Math 

Variable 

Code 

Variable 

Default 

Value 

Mean  RDE  throughflow  diameter 

Dm 

D_m 

140  mm 

RDE  thickness 

T 

tau 

10  mm 

RDE  height 

hRDE 

Ht 

177  mm 

Total  pressure  entering  RDE 

P t3  ,air 

Pt3_air 

58.784  psia 

Pt3,fuel 

Pt3_fuel 

58.784  psia 

Total  temperature  entering  RDE 

Pt3,air 

Tt3_air 

540  R 

Pt3,fuel 

Tt3_fuel 

540  R 

Oblique  shock  angle  relative  to  flow 

@ shock 

theta_shock 

60° 

Number  of  detonation  waves 

'H'waves 

n^wave 

1 

Deflagration  flame  speed  at  standard  conditions 

Vdefl,std 

V_f lame 

10  ft/s 

Injector-to-throat  area  ratio 

■A inj,air 

As 

AinjqA3_air 

0.1801 

s^inj  ,fuel 

As 

AinjqA3_fuel 

0.0199 

Station  3  flow  coefficient 

C'W3,fwd,air 

cw_fwd_air 

1 

C'W3,bwd,air 

cw_back_air 

0 

C'W3,fwd,fuel 

cw_fwd_fuel 

1 

C'W3,bwd,fuel 

cw_back_fuel 

0 

Ambient  pressure 

Po 

PO 

14.7  psia 

Combustion  chamber  area  contraction 

^8 

A8qA3 

1 

RDE  exit  flow  coefficient 

C\V8 

Cw8 

1 

Nozzle  area  ratio 

Ag 

As 

A9qA8 

1 

Nozzle  stream  thrust  coefficient 

Cs 

Cs 

1 

Axial  Mach  number  at  throat,  path  a 

Masy 

M8y_det 

1 

Axial  Mach  number  at  throat,  path  b 

MbSy 

M8y_det_shk 

1 

Axial  Mach  number  at  throat,  path  d 

Md8y 

M8ynnix 

1 

Detonation  lateral  area  relief 

(feL 

A2qAl_C J 

1 
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Initial  Calculations 


Before  the  NPSS  model  implements  the  Kaemming  model  illustrated  in  Fig.  11,  several 
parameters  are  calculated  based  on  the  input  conditions  from  Table  6.  In  the  Kaemming 
model,  these  parameters  are  automatically  calculated  (based  on  the  input  parameters)  be¬ 
cause  Microsoft  Excel  automatically  updates  linked  cells;  however,  NPSS  does  not  automat¬ 
ically  update  linked  cells,  so  the  values  must  be  manually  updated. 

Geometric  Parameters. 

The  first  group  of  calculations  are  for  geometric  parameters.  The  inner  and  outer  di¬ 
ameters  are  calculated  using  Eq.  11  and  Eq.  12,  respectively.  The  effective  circumference 
per  detonation  wave  is  calculated  using  the  mean  diameter  from  Eq.  13.  The  cross-section 
through-flow  area  per  detonation  wave  is  calculated  using  Eq.  14.  The  air  and  fuel  injector 
areas  are  calculated  by  multiplying  A3  by  the  ratios  A"in^ir  and  Atn^fuel ;  respectively,  from 
Table  6.  In  the  default  case,  nwave  =  1. 


Di  —  Dm  T 

(ip 

D0  —  Dm  T  t 

(12) 

1 vDm 

^e//  — 

J  J  77 

1  Lwave 

(13) 

\  '  '"wave  / 

(14) 

Gas  Properties  for  Input  to  NPSS. 

The  second  group  of  calculations  are  for  gas  properties.  In  the  NPSS  model,  gas  properties 
are  calculated  based  on  two  variables:  fuel  and  gas_prop_f lag.  The  fuel  parameter  indicates 
the  fuel  used,  either  hydrogen  H2  or  ethylene  C2H4.  Table  7  summarizes  the  gas  properties  for 
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the  air,  hydrogen,  and  ethylene.  Gas  properties  for  the  fuel  and  air  reactants  as  well  as  the  post¬ 
detonation  mixed  products  are  calculated  ahead  of  time  outside  of  the  NPSS  model  using  either 
known  values  or  the  NASA  computer  program  Chemical  Equilibrium  with  Applications  (CEA). 

The  variable  gas_prop_f  lag  determines  if  the  gas  properties  for  the  reactants  are  taken  from  a 
lookup  table  or  calculated  [4] .  Kaemming  used  the  tabulated  lookup  values  summarized  in  Table  8 
for  hydrogen  from  CFD  simulations  [3;  14]  and  for  ethylene  from  CEA  [4;  15]. 


Table  7.  Gas  properties  for  air,  hydrogen,  and  ethylene 


Parameter 

Air 

Hydrogen 

Ethylene 

Molecular  weight,  MW  ( kg/kmol ) 

28.965 

2.016 

28.05 

Ratio  of  specific  heats,  7 

1.400 

1.405 

1.2313 

Specific  gas  constant,  R  (  kgJ_K  ) 

287.052 

4124.236 

296.416 

Stoichiometric  fuel-to-air  ratio,  FARstoich 

- 

0.0292 

0.0676 

Heat  of  combustion,  AH 

- 

3550 

3181 

Table  8.  Tabulated  gas  properties  for  reactants 


Parameter 

Hydrogen 

Ethylene 

Ratio  of  specific  heats,  7 r 

1.4256 

1.3690 

Specific  gas  constant,  Rr  (^kqLK'j 

397.500 

288.807 

Kaemming  calculates  the  reactant  properties  using  the  inlet  conditions  from  Table  6  to  find 
the  fuel-to-air  ratio  at  the  detonation  and  then  mole-averages  the  gas  properties  from  Table  7  [4]. 
Kaemming  assumed  the  injectors  are  choked,  so  the  mass  flow  parameter  for  the  air  and  fuel  are 
calculated  using  Eq.  15  and  Eq.  16,  respectively  [7]. 


MFPnjr  = 


' lairdc  I  .  'Yair  1 

1 


'Yair  +  1 


^(l  'Yair ) 


MFPfuei  =  „  /7£Mdgc  ( 1  +  7fUel  1 


R 


fuel 


'Yfuel  +  l 
K'-Vuel) 


(15) 


(16) 
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The  mass  flow  parameter  is  related  to  the  mass  flow  by  Eq.  17  [7:  Eq.  1.3].  The  detonation  fuel- 
to-air  ratio  was  calculated  as  a  ratio  of  the  fuel  mass  flow  rate  to  the  air  mass  flow  rate  using 


Eq.  18. 


/ Ainj,fuel\  ( MFPfUel\  I  Pt3,air 

\  Ainj,air  )  \  M F Pair  )  y  Tt:i  Juei 
The  mole  fraction  for  a  species  i  is  defined  using  Eq.  19  [18]. 

m.j 

v.  -  MWi 
^  i  m-i 

Mi  MWi 


F  AR.rift  — 


mfuel 


m„ 


PfA.f  at  l 

Pt3,air 


(17) 

(18) 


(19) 


The  fuel-to-air  ratio  is  used  to  calculate  the  mole  fraction  for  air  and  fuel  using  Eq.  20  and  Eq.  21, 
respectively. 


Xntr  — 


MW n 


,  MWa 


+ 


FARdet 

MWfuel 


(20) 


(_L+W)  <21> 

v  Mel/  \  MWair  -1-  MWfuei  ) 

Kaemming  uses  the  mole  fractions  to  take  a  weighted  average  of  the  molecular  weight  and  ratio  of 
specific  heats  using  Eq.  22  and  Eq.  23,  respectively  [4;  19]. 


MWr  =  MWairXair  +  MW fuelX  fuel 


(22) 


TR  —  1 air  X  air  +  'J fuelX  fuel 


(23) 


The  molecular  weight  is  used  to  calculate  the  specific  gas  constant  using  Eq.  24. 


Rr 


Ru 

MWr 


(24) 


In  both  the  tabular  lookup  (Table  7  and  Table  8)  and  calculated  cases  (Eqs. 
properties  for  the  products  are  taken  from  a  lookup  table  based  on  CEA  results  [4; 


15-24),  the  gas 
15].  Kaemming 
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uses  fixed  values  for  the  products  because  the  properties  cannot  be  calculated  a  priori  [4;  15]. 

The  assumption  of  constant  Rp  and  7 p  is  a  valid  approximation  to  actual  properties,  even  with 
variations  in  pressure  and  temperature  [2],  Table  9  summarizes  the  values  used  for  the  products. 


Table  9.  Tabulated  gas  properties  for  products 


Parameter 

Hydrogen 

Ethylene 

Ratio  of  specific  heats,  7 > 

1.2412 

1.2300 

Specific  gas  constant,  Rp 

347.700 

298.191 

Entrance  Flow  Model 


In  the  NPSS  model,  the  subroutine  iter8  initializes  the  iteration  for  the  pressure  at  station 
3.2  and  ensures  the  fuel-to-air  ratios  are  appropriately  updated.  Both  the  Kaenuning  and  NPSS 
models  assume  separate  injectors  for  the  air  and  fuel  streams;  the  parameters  Am^™r  and  Aznrfuel 
from  Table  6  describe  the  proportion  of  the  total  area  made  up  of  air  or  fuel  injectors,  respectively 
[4],  Additionally,  both  models  allow  for  forward  (injection)  and  backward  (blowback)  flow  for  the 
air  or  fuel  injectors  [4];  the  station  three  flow  coefficients  (C\y 3)  from  Table  6  govern  the  forward 
or  backward  direction  of  flow  for  the  air  and  fuel.  In  the  default  setup,  all  flow  is  assumed  to  be 
forward  flow  (Cw;$jwc[  =  1)  with  no  backflow  (CW2,bwd  =  0)  [4], 

Kaemming  finds  a  minimum  value  of  P3.2  based  on  isentropic  expansion  from  the  injector  into 
the  RDE  chamber  (i.e.,  station  3*  to  substation  3.2)  [4],  The  amount  of  expansion  is  based  on  the 
total  area  of  injectors  and  the  total  area  in  the  RDE  chamber,  given  by  Eq.  25.  The  expansion  is 
multiplied  by  the  mass-averaged  forward  mass  flow  coefficient,  given  by  Eq.  26. 


A-t 


Ainj, total  A^r,air  +  Anpfuel 


(25) 


A3 


Cw3Jwd,eff  = 


C'W‘i)fwd,air  T  F  ARdet Cw3,fwd,fuel 


1  +  FAR, 


det 


(26) 


The  maximum  Mach  number  associated  with  the  area  expansion  A  3 —  is  found  using  the  function 


inj,  total 
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Table  10.  Grid  convergence  study  for  discretized  entrance  flow  model  with  percent  change 
based  on  100-step  solution 


Steps 

Time  (sec) 

Specific  Impulse  (sec) 

%  Change 

50 

3 

2922.28 

1.6% 

100 

4 

2970.43 

- 

200 

6.5 

2986.29 

-0.53% 

Msup_A.  The  Mach  number  is  used  to  find  the  minimum  pressure  at  station  3.2  using  the  function 
PtqP  and  the  total  pressure  Ptz  entering  the  RDE  from  Table  6. 

The  entrance  flow  model  calculates  the  thermodynamic  properties  of  the  flow  entering  the  RDE, 
as  seen  in  Fig.  11  [4].  In  the  NPSS  element,  the  entrance  flow  model  calculations  are  performed  in 
the  f  ill_calc  subroutine.  The  duration  of  each  cycle  is  calculated  using  Eq.  27  [4].  The  calculation 
of  the  detonation  velocity  V)let  is  described  later. 


t 


cycle 


Ceff 

Vdet 


(27) 


Each  cycle  is  discretized  into  100  steps.  Table  10  summarizes  a  grid  convergence  study  showing 
the  number  of  discretized  steps,  time  to  complete  the  entire  solution,  calculated  specific  impulse, 
and  percent  change  from  the  100-step  solution.  The  specific  impulse  calculation  used  the  default 
input  values  from  Table  6.  In  all  three  cases  summarized  in  Table  10,  the  specific  impulse  does 
not  significantly  change,  suggesting  that  the  solution  has  reached  a  stable  solution.  The  use  of  100 
steps  provides  a  balance  of  solution  precision  and  computational  time. 

The  f  ill_calc  subroutine  uses  P3.2  from  the  iter8  subroutine  in  the  empirical  plenum  model 
from  Eq.  1  and  Eq.  2  to  model  the  pressure  profile  inside  the  combustion  chamber  [4].  The  modeled 
pressure  determines  if  the  flow  is  moving  forwards  or  backwards  by  comparing  the  modeled  pressure 
to  the  pressure  of  the  air  and  fuel  from  Table  6  [4].  The  direction  of  the  flow  dictates  an  increase 
or  decrease  in  air  or  fuel,  and  the  amount  of  air  or  fuel  mass  flow  is  calculated  using  Eq.  28  with 
the  appropriate  values  based  on  Table  11.  The  mass  flow  rate  is  multiplied  by  the  length  of  the 
discretized  time  step  to  calculate  the  mass  of  the  air  and  fuel  used.  The  injection  Mach  number 
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Minj  has  an  upper  limit  of  one  for  the  choked  condition  [4].  The  properties  at  substation  3.4  are 
calculated  outside  of  the  f  ill_calc  subroutine  and  is  described  later. 


m = CwPAM 


(28) 


Table  11.  Parameters  used  in  air  and  fuel  All  calculations  based  direction  of  flow  using  values 
from  Table  6 


Forward 

Backward 

Air 

C'W  C'W,fwd,air 

M  =  f{Pt3, air/P3A ) 

P  =  f{Pt3,air,  M) 

Tt.air  Tt3,air 

T  =  f(Tt3,air,  M) 

A  —  A.  .  . 

■  1  -  *-inj.air 

7  =  "lair 

MW  =  MWair 

R  =  ^ u 

MW„.ir 

C’W  Cw,bwd,air 

M  =  f{P3A/Pt3,air ) 

P  =  f(P3,  M) 

Tt,air  =  Tt,a3A 

T  =  f(Tta3A,  M) 

A  —  A.  .  . 

-  1  -  '-i.nj.air 

1  =  1P 

MW  =  MWP 

R  =  Ru 

1  MWp 

Fuel 

Cw  =  CwjwcLJuel 

M  =  f(Pt3Juel/P3A) 
P  =  f(Pt3Juel,  M) 

Ttjuel  Tt3,fuel 

T  =  f(Tt3Juel,  M) 

A  Ainjjuei 

T  'y  fuel 

MW  =  MW fuel 

R  =  ^ u 

MWfnel 

C'W  =  Cw,bwd,fuel 

M  =  f(P3A/Pt3,fuel) 
P  =  f{P3r  M) 

Tt.fuel  Tt,a3A 

T  =  f(TUa3. 4,  M) 

A  Ainjjuei 

7  =  7  p 

MW  =  MWp 

R  =  Ru 

1L  MWp 

The  fill.calc  subroutine  tracks  masses  of  the  air  and  fuel  throughout  each  iteration.  The 
tracked  masses  of  air  and  fuel  are  used  to  calculate  the  total  temperature,  ratio  of  specific  heats, 
and  molecular  weight  of  the  mixture  using  mass- averaging  based  on  Eqs.  29-  Eq.  31,  respectively, 
using  values  from  Table  6  and  Table  11  [4;  14]. 

m  AlairTt  dir  T  AlfUelPt ,fUel  /nn\ 

Tt,  3.2  =  - r - —r - - —  (29) 

'W'air  Wlfuel 
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73.2  = 


MW’s.  2  = 


MWair)  laiT  +  T/uei 


+ 


T'fuel 


MWair  ^ 


MWa 


~  )  fuel 


+ 


fuel 


MWair  MWfuel 

The  static  pressure  mass  flow  parameter  definition  is  given  in  Eq.  32  [7:  Eq.  1.4]. 


(30) 


(31) 


MFP=7^§  =  Wf  (X  +  ^M2 


(32) 


Squaring  Eq.  32  and  rearranging  the  result  leads  to  Eq.  33 


7“  1  fn,2\2 


(mT  +  m2+(A)2(^, 


19c 


=  0 


(33) 


The  quadratic  equation  is  used  to  find  the  the  Mach  number  given  in  Eq.  34. 


M  = 


\ 


-2 

7-1 


+4fAY7  ^ 

V73c(7-1) 


(34) 


Kaemming  uses  Eq.  34  to  calculate  the  Mach  number  at  station  3.2,  as  seen  in  Eq.  35. 


M:s/2  = 


\ 


-2 


73.2  -  1 


+ 


\ 


73.2  -  1 


A  ( 'W'dir  rilfuel  \ 

+  4  pUs  ) 


MW3.2  )  ^*’3.2 


73.2^(73.2  -  1) 


(35) 


The  total  temperature  and  Mach  number  at  station  3.2  are  used  to  calculate  the  static  temper¬ 
ature.  The  station  3.2  properties  are  used  to  calculate  the  velocity  of  the  flow  using  Eq.  36  based 
on  the  definition  of  mass  flow  rate. 


^3.2 


( rilair  T  TTlfuel) 


(  A/U'x..  \ 

\  T3.2R.U  ) 


^3 


(36) 


Kaemming  assumes  the  deflagration  flame  speed  is  proportional  to  density  [20]  based  on  Eq.  37. 
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The  default  flame  speed  of  10  ft/ sec  from  Table  6  is  at  a  standard  pressure  of  Pst,d  =  14.7  psia 
and  a  standard  temperature  of  Tstd  =  518.7  R  [4]. 


T/  lame,  local  —  V flame, std 


P 


P. 


std 


(37) 


Kaemming  tracks  the  deflagrated  mass  by  comparing  the  local  flame  speed  to  V3.2.  If  V3.2  is  less 
than  the  local  flame  speed,  all  of  the  air  and  mass  flow  at  the  given  time  step  is  deflagrated  [4],  If 
V3.2  is  greater  than  the  local  flame  speed,  only  a  proportion  of  the  air  and  mass  flow  is  deflagrated 
based  on  Eq.  38  [4], 


A mdefi  =  (rhair  +  riifuei)  A t 


V3, 


V f  lame,  local 


(38) 


Additionally,  Kaemming  uses  V3.2  to  calculate  the  fill  height  of  the  flow  during  the  cycle  using 
Eq.  39  [4], 

A  hfni  =  V3.2A  t  (39) 


At  each  of  the  100  steps,  the  fill.calc  subroutine  calculates  and  stores  the  following  param¬ 
eters  in  arrays:  mass  of  air  used,  mass  of  fuel  used,  and  fill  height  [4],  All  three  parameter  arrays 
correspond  to  values  in  an  array  for  the  total  mass,  so  the  three  parameter  arrays  can  be  seen  as 
functions  of  the  total  mass.  The  fill_calc  subroutine  calculates  the  total  mass  of  deflagrated 
products  by  summing  the  incremental  calculations  from  Eq.  38  and  uses  the  function  OneDNtrp  to 
interpolate  along  the  total  mass  array  and  find  corresponding  values  for  the  mass  of  air  deflagrated, 
mass  of  fuel  deflagrated,  and  deflagration  fill  height  [4].  The  percentage  of  the  flow  composed  of 
deflagration  is  calculated  by  dividing  the  deflagration  mass  by  the  total  mass  at  the  end  of  the  100 
steps,  as  seen  in  Eq.  40  [4].  The  fuel-to-air  ratio  for  the  deflagrated  flow  is  calculated  using  Eq.  41. 


ry  (  fuel, defl  T  fTl air,defl\ 

A defl  =  — - - - — 

V  mtot  J 


(40) 


FARdefl  = 


m  fuel, defl 


(41) 


YYl air, defl 

Kaemming  assumes  the  remainder  of  flow  not  deflagrated  is  detonated  [4].  The  detonation 
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height  is  calculated  by  subtracting  the  deflagration  height  (previously  found  by  interpolation)  from 
the  total  fill  height  at  the  end  of  the  100  steps.  The  detonation  height  is  divided  by  the  height 
of  the  RDE  to  calculate  the  parameters  from  the  empirical  models  from  Fig.  9  and  Fig.  10.  The 
percent  of  the  flow  composed  of  detonation  is  calculated  using  Eq.  42.  The  detonation  fuel-to-air 
ratio  originally  estimated  in  Eq.  18  is  updated  using  Eq.  43. 

Zdet  =  1  —  Zdefl  (42) 

FARdet  =  mfuel’tot  ~  (43) 

Wlair,tot  ^ air,defl 

Kaemming  calculates  the  average  mass  flow  rate  using  Eq.  44.  The  average  mass  flow  rate 
calculated  in  Eq.  44  is  “w3”  in  Fig.  11. 


ril3,avg 


mtat 


t  cycle 


(44) 


In  the  Kaemming  model,  the  fill  cycle  calculations  are  used  to  update  the  air  valve/injector 
model,  as  seen  in  Fig.  11  [4].  In  the  NPSS  element,  the  subroutine  Inject.calcs  performs  the 
the  air  valve/injector  model  calculations.  The  fill  cycle  calculations  are  based  on  post-combustion 
properties  while  the  air  valve/injector  model  calculations  are  based  on  pre-combustion  properties. 

The  injector  Mach  number  Mmj  for  the  air  and  fuel  are  calculated  based  on  the  ratio  of  their 
respective  total  pressures  from  Table  6  to  the  T3.2  for  the  current  iteration.  The  air  and  fuel 
flow  Mach  numbers  have  an  upper  limit  of  one  for  choked  flow.  The  Mach  numbers  are  used  to 
calculate  the  pressure  Pinj  and  temperature  Tjnj  of  the  air  and  fuel  flows  inside  the  injector.  The 
Mach  number,  pressure,  and  temperature  are  used  to  calculate  the  mass  flow  for  the  air  and  fuel 
using  Eq.  45  with  the  parameters  from  Table  12  [4].  Table  6  summarizes  the  default  values  for  the 
parameters  in  Table  12. 

minj  =  CwPAmJ H  (45) 
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Table  12.  Parameters  used  for  mass  flow  calculations 


Parameter 

Air 

Fuel 

Cw 

Cw3Jwd,air 

C'W3,fwd,fuel 

P 

P 

1  in  j, air 

Pin  j,  fuel 

A 

A 

■^inj^air 

-A-inj,fuel 

M 

M-inj,air 

MinjJUel 

7 

air 

'y fuel 

R 

Pair 

Pfuel 

T 

T 

-1-  mj,air 

Pin  j,  fuel 

The  pressure  and  temperature  streams  are  combined  to  get  the  mass-averaged  total  pressure 
and  total  temperature  using  Eq.  46  and  Eq.  47,  respectively  [4]. 


Pi 


tS,sum  — 


P t3, air  P ARdetPt3,fuel 


1  +  FAR, 


=  P 


det 


t3 


(46) 


Pt3,sum  — 


C P:airPt3,  air  T  F AR,(ietC pjuel Ft 3 , / u e l 


(1  +  F  ARdet)Cp,R 


=  Tt  3 


(47) 


Kaemming  calculates  the  Mach  number  at  station  3.2  based  on  the  steady-state  mass  flow  rate, 
static  pressure,  and  total  temperature  using  Eq.  48  [4],  The  calculation  for  M3.2  in  Eq.  35  is  for  a 
variable  used  internally  in  the  fill.calcs  subroutine. 


M3. 2 


1 

2 


2  i  f  ^  \  1  1  Flinj,air  T  Pinj,fuel  \  f  2 RRTtZ,sum  \ 

1r  ~  1  V  v /'/?  ~  V  V  P3.2A3  J  -  1 )/ 


(48) 


The  Mach  number  and  isentropic  relations  functions  are  used  to  calculate  the  station  3.2  total 
pressure  Pt, 3.2  and  static  temperature  T3.2  from  P3.2  and  7)3, .sum ,  respectively.  The  axial  velocity 
at  station  3.2  is  calculated  using  Eq.  49  [4], 


V3.2 y  =  M3. 2  yj ''1rRr9cTz.2 


(49) 
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Thermodynamics  Model 


The  exit  flow  model  calculates  the  thermodynamic  properties  of  the  flow  exiting  the  RDE,  as 
seen  in  Fig.  11  [4].  In  the  NPSS  element,  the  exit  flow  model  calculations  are  performed  in  the 
Thermo.parametric  subroutine.  The  exiting  flow  is  composed  of  the  four  paths  a-d  in  Fig.  7. 
Kaemming  calculates  the  thermodynamic  properties  for  the  flow  along  each  pathline  to  the  exit  of 
the  RDE  and  then  expands  each  stream  to  the  exhaust  based  on  throat  Mach  numbers  defined  in 
Table  6  [4],  After  the  properties  for  the  separate  streams  are  calculated  all  the  way  to  ambient 
conditions,  properties  such  as  total  temperature,  specific  thrust,  and  specific  impulse  are  mass- 
averaged  [4], 


Path  a:  Detonation  Only. 

The  post-detonation  (station  a3.4  in  Fig.  7)  properties  are  calculated  using  properties  for  the 
products  and  enthalpy  addition.  The  total  temperature  at  substation  a3.4  is  calculated  by  assuming 
enthalpy  addition  to  the  station  3  properties  based  on  Eq.  50  [4]. 

Tt,a3A  =  (Cp,RTt3,sum  +  l  +  FARdet )  (5°) 

The  non-dimensional  heat  addition  term  is  calculated  using  Eq.  51  [4]  (see  Appendix  A,  Eqs.  A. 21- 
A. 24).  Substation  3.2  properties  CpF  T3.2  are  used  as  the  pre-detonation  properties  Cp,\  T\  from 
Eq.  A. 24. 

AH  (1  +  F ARdet)  \ 

“ =  (51) 

The  detonation  Mach  number  M(jet  is  calculated  using  the  q  term  from  Eq.  51  in  the  function 
CJ_mach_area.  The  detonation  velocity  is  calculated  using  Eq.  52. 


Vdet.  —  Mdet\/  1r9cRrTz.2 


(52) 


The  pressure  and  temperature  change  across  the  detonation  with  lateral  area  relief  are  calculated 
using  Eq.  53  and  Eq.  54,  respectively.  By  default,  no  lateral  area  relief  is  assumed  (i.e.,  = 

1);  further  research  is  required  to  find  appropriate  values  for 


Pa'S  A  =  P'S. 2 


7 RMdet  +  1 

7p('feh.  +  1 


CJ 


(53) 


TaSA  ~  T3.2 


Cp,R  \  f  1  +  +  1 


CPA 


7P  +  1 
2 


(54) 


The  static  temperature  and  static  pressure  are  used  to  calculate  the  entropy  change  to  station  3.4 
based  on  Eq.  55. 


^^3.4  =  [Cp,p  ln(Ta3.4)  —  CP,Rln(T3.2)] 


(55) 


Kaemming  changes  the  frame  of  reference  from  the  laboratory  frame-of-reference  to  the  det¬ 
onation  wave  frame-of-reference  using  the  velocity  triangles  in  Fig.  5  [4],  In  the  detonation  wave 
frame-of-reference,  the  ZND  detonation  model  is  valid  [4;  5;  11].  The  velocity  of  the  flow  coming 
into  the  detonation  wave  in  the  detonation  frame-of-reference  is  calculated  as  a  vector  sum  of  the 
axial  velocity  and  detonation  wave  velocity,  as  seen  in  Eq.  56. 


(56) 


Kaemming  assumes  the  circumferential  Mach  number  after  the  detonation  in  the  detonation  wave 
frame-of-reference  Msax,w  is  one  [4] .  The  circumferential  velocity  in  the  detonation  wave  frame-of- 
reference  is  calculated  using  Eq.  57 


T 3Ax,w  —  P^3Ax,w  \! 7  P fir P /  ’ Tj . 4 


(57) 


Kaemming  converts  the  static  pressure  and  static  temperature  calculated  in  Eq.  53  and  Eq.  54 
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to  stagnation  in  the  detonation  wave  frame-of-reference  (which  is  equivalent  to  accelerating  the  flow 
to  the  detonation  velocity  in  the  laboratory  frame-of-reference)  to  calculate  the  total  temperature 
in  the  laboratory  frame-of-reference  using  Eq.  58  [4;  11].  The  total  temperature  calculated  in 

Eq.  58  was  compared  to  the  total  temperature  calculated  in  Eq.  50  to  ensure  agreement  [4],  The 
total  temperature  is  used  to  calculate  the  total  pressure  using  isentropic  relations,  as  seen  in  Eq.  59 


[21]- 


Tt,a3A  ~  Ta 3.4  + 


V2 

v 3Ax.w 


-vL  +  v* 


3.2  y 


2  Cpj 


Rt,a3A  —  Pa3A 


/  Tttg3A 
V  T<x3A 


.  'if p 

7  p  — 1 


(58) 

(59) 


The  total  temperature  and  total  pressure  calculated  from  Eq.  58  and  Eq.  59  are  used  to  calculate 
the  total  entropy  change  based  on  the  total  properties  using  Eq.  60. 


ASt,a3.4  =  [Cp,p  In  (Ttta3A)  -  CptR  ln(T3.2)] 


Rp  In 


f  Pt,a3A 


(60) 


The  flow  traveling  along  path  a  continues  to  the  exit  of  the  RDE  chamber  with  substation 
3.4  properties,  as  seen  in  Fig.  7.  After  exiting  the  the  RDE,  the  flow  is  expanded  to  the  throat 
(station  8),  exit  (station  9),  and  ambient  using  isentropic  relations.  The  circumferential  velocity 
calculated  from  the  empirical  model  described  in  Eq.  6  is  assumed  to  remain  constant  throughout 
the  expansion  process  because  no  energy  is  extracted  from  the  rotating  flow  [4;  7].  The  Mach 
number  at  the  throat  is  calculated  based  on  the  user-defined  axial  throat  Mach  number  (assumed 
to  be  choked  by  default  from  Table  6)  and  the  empirically  modeled  circumferential  velocity,  as  seen 
in  Eq.  61  [4], 


Ma8 


MaSy 


(61) 


(y  1  —  (Circumferential  Velocity  Detonated  Flow)2 

The  Mach  number  at  the  throat  is  used  to  calculate  the  pressure,  temperature,  and  velocity  of 
the  flow  at  the  throat  using  the  isentropic  relations  functions  [21].  Kaemming  uses  the  user-defined 
nozzle  stream  correction  factor  (Cs)  and  empirical  correction  factors  from  Fig.  10  to  calculate  the 
specific  thrust  using  Eq.  62  [3;  4;  14].  The  specific  thrust  is  based  on  the  stream  thrust  function 
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[22:  Eq.  2.92]  divided  by  the  mass  flow  rate  equation  from  Eq.  28.  Kaemming  uses  the  mass 
correction  term  (Eq.  9)  and  momentum  correction  term  (Eq.  10)  to  correct  the  mass  flow  rate  and 
specific  thrust,  respectively  [15]. 


In  rocket  engines,  the  specific  impulse  is  given  as  the  ratio  of  the  thrust  to  the  weight  flow  of 
fuel  and  oxidizer  [22],  In  air-breathing  engines,  the  specific  impulse  is  given  as  the  ratio  of  the 
thrust  to  the  weight  flow  of  the  fuel  alone  [22],  In  each  case,  the  specific  impulse  is  based  on 
what  the  vehicle  is  required  to  carry  (i.e.,  both  the  fuel  and  oxidizer  in  the  rocket  or  only  fuel  in 
the  air-breathing  engine).  The  specific  impulse  calculated  in  this  thesis  is  the  specific  impulse  for 
air-breathing  vehicles.  The  specific  impulse  for  rockets  can  be  found  by  multiplying  the  specific 
thrust  by  ■  Note  that  the  specific  impulse  for  rockets  will  have  the  same  magnitude  as  the 
specific  thrust.  The  specific  impulse  for  air-breathing  engines  is  calculated  using  Eq.  63  [3;  4;  14]. 


_  (F\  fl  +  FARdet\  fgA 
'a8  U/asV  FARdet  ){g0J 


(63) 


The  flow  at  the  throat  is  expanded  to  the  exit  (station  9).  The  axial  Mach  number  at  the  exit 
Ma 9y  is  calculated  using  the  function  expans ionunode.  The  total  temperature  between  stations  8 
and  9  is  conserved  [7;  22],  so  it  is  used  to  calculate  a  parameter  analogous  to  the  circumferential 
Mach  number,  as  seen  in  Eq.  64. 


n2  v&x 

at )  lfP9cRpTt,a3A 


(64) 


The  station  9  Mach  number  is  calculated  using  Eq.  65,  which  corrects  for  the  total  temperature 
used  in  Eq.  64. 


Ma9  — 


(65) 


The  exit  Mach  number  Mag  and  isentropic  relations  functions  are  used  to  calculate  the  temperature, 
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pressure,  and  speed  of  sound  at  station  9  [4;  21].  The  station  9  velocity  is  calculated  as  the  product 
of  the  station  9  Mach  number  and  speed  of  sound.  The  axial  velocity  component  of  the  exit  flow 


is  calculated  using  Eq.  66  [4], 


V«*y  =  -  V?9x 


(66) 


The  specific  thrust  and  specific  impulse  are  calculated  using  Eq.  67  and  Eq.  68  [3;  4;  14]. 


F\  Cs  &)  (l  +  VM&,)  ~  % 

’hJ‘°  (aMdgfe) 


T  _  (  F\  1  +  F ARdet  \  (  gc 

SP’a9~Vn)aA  FARdet  ) 


(67) 


(68) 


a9  \  ±  -^±Ldet  /  \90, 

To  expand  the  flow  to  the  ambient  conditions,  the  Mach  number  M0  amb  is  calculated  using 
the  pressure  ratio  Ptp^ 4  and  the  function  M_PtqP.  The  Mach  number  is  used  to  calculate  the 
temperature,  area  ratio,  speed  of  sound,  and  velocity  using  isentropic  relations  functions  [4;  21]. 
The  constant  circumferential  velocity  is  again  used  to  remove  the  circumferential  component  from 
the  total  velocity  similar  to  Eq.  66,  resulting  in  the  axial  velocity  component  Vay>amb-  The  ideal 
specific  thrust  is  calculated  using  Eq.  69  [4;  22],  The  associated  specific  impulse  is  calculated  using 
Eq.  70. 

/  r\  \  ■ .  , 

(69) 
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ay,amb 


*■ sp,amb 


F 


rn 


a,amb 


a,amb 


9c 


1  +  FARdet  \  (  9c 


FARdet  J  V  9o 


(70) 


Path  b:  Detonation  and  Shock. 


In  Fig.  7,  the  flow  traveling  along  path  b  passes  through  the  same  detonation  between  substation 
3.2  to  63.4  as  the  flow  traveling  along  path  a  from  substation  3.2  to  a3.4.  Therefore,  the  properties 
at  substation  63.4  equal  those  at  substation  a3.4.  Unlike  the  flow  traveling  along  path  a,  the  flow 
traveling  along  path  6  passes  through  an  oblique  shock  between  substation  63.4  and  substation  63.6. 
Kaemming  calculates  the  Mach  number  entering  the  detonation  wave  in  the  detonation  wave 
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frame-of-reference  using  Eq.  71  [4]. 


Mdet,w 


Vdet 

V  iPfjcMpT'i.A 


(71) 


The  normal  component  of  M^etw  is  calculated  using  the  user-defined  angle  of  the  shock  relative 
to  the  flow  Oghock  from  Table  6.  The  shock  angle  is  used  to  determine  the  normal  component  of 
the  Mach  number.  The  normal  component  of  the  Mach  number  is  used  with  the  normal  shock 
functions  to  find  the  temperature  and  pressure  at  station  63.6  [21].  The  greater  the  shock  angle, 
the  greater  the  normal  component  of  the  Mach  number,  and  the  greater  the  change  in  temperature 
and  pressure  across  the  oblique  shock  [21].  The  entropy  rise  across  the  shock  is  calculated  using 
Eq.  72  [4],  The  entropy  rise  across  the  shock  calculated  in  Eq.  72  is  added  to  the  entropy  rise  across 
the  detonation  calculated  in  Eq.  55  to  find  the  total  entropy  rise  along  path  6  [4], 

A5b3.6  =  Cp.p  In  -  RP  In  (72) 


The  flow  traveling  along  path  b  continues  to  the  exit  of  the  RDE  chamber  with  substation 
63.6  properties,  as  seen  in  Fig.  7.  After  exiting  the  the  RDE,  the  flow  is  expanded  to  the  throat 
(station  8),  exit  (station  9),  and  ambient  using  isentropic  relations  (see  expansion  of  path  a  flow 
for  process)  [4]. 


Path  c:  Deflagration. 

Kaemming  assumes  the  ratio  of  specific  heats  and  gas  constant  are  linear  functions  of  the  fuel- 
to-air  ratio  [4].  If  the  detonation  fuel-to-air  ratio  is  greater  than  the  deflagration  fuel-to-air  ratio, 
the  ratio  of  specific  heats  and  gas  constant  at  station  c3.4  are  calculated  using  Eq.  73  and  Eq.  74, 
respectively  [4],  If  the  deflagration  fuel-to-air  ratio  is  greater  than  the  detonation  fuel-to-air  ratio, 
the  ratio  of  specific  heats  and  gas  constant  at  station  c3.4  are  calculated  using  Eq.  75  and  Eq.  76, 
respectively  [4], 

7c3.4  =  7 air  +  (7 P  ~  lair)  (73) 
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PC3.4  =  Rair  +  (Rp  ~  Rair )  f  pAR  ^ 


7c3.4  =  7  P 


Rc3A  =  Rp 


(74) 

(75) 

(76) 


The  ratio  of  specific  heats  and  gas  constant  are  used  to  calculate  molecular  weight  and  specific  heat 
capacity  using  Eq.  77  and  Eq.  78,  respectively  [4;  21]. 


MW*.  4  = 


Ru 

Rc3A 


Cp,c3A  ~ 


7c3.4 


Rc3A 


(77) 


(78) 


,7c3.4  -  1, 

Kaemming  calculates  the  total  temperature  at  station  c3.4  using  Eq.  79,  where  FARmin  is  the 
the  lower  of  FAR^efi  and  F AR^et  [4] . 


Tt,c3A  — 
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Cp,c3.4 


Cp,RTt^sum  + 


AH 


1  +  . 


(FARn 


FARdet )  \  FARdet 


(79) 


Kaemming  assumes  the  deflagration  is  a  constant  pressure  process,  so  the  pressure  at  substation 
c3.4  is  the  same  as  the  pressure  at  substation  3.2  [4],  The  total  temperature  T))C3.4  is  used  with  the 
axial  velocity  V3.2 y  from  Eq.  36  to  calculate  the  axial  Mach  number  at  station  c3.4  using  Eq.  80. 


Mc3Ay  — 


V2 
*4.2  y 


\|  Ac3A9cRc3ATtc3A  ~  (^f1)  V322y 


(80) 


The  Mach  number  is  then  used  to  calculate  the  static  temperature  using  isentropic  relations  [21]. 
The  entropy  rise  after  deflagration  is  calculated  using  Eq.  81  [4]. 


A5C3.4  =  [Cp)C3.4  ln(Tc3.4)  —  C'piijln(T3.2)]  - 


Rc3A  la  (  -  Rr  hi 


Pn 


Pn 


(81) 


The  flow  traveling  along  path  c  passes  through  an  oblique  shock  between  substations  c3.4  and 
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c3.6,  as  seen  in  Fig.  7.  Kaemming  changes  the  frame  of  reference  from  the  laboratory  to  detonation 
wave  frame-of-reference  based  on  Fig.  5  in  order  to  use  normal  shock  relations  [4;  21].  The  total 
velocity  of  the  flow  is  given  by  Eq.  56  and  the  circumferential  component  is  given  by  from 
Eq.  52  [4].  The  component  of  the  Mach  number  normal  to  the  oblique  shock  wave  is  calculated 
using  Eq.  82  [4], 


Ti c3.4x  — 


Ec3. 


Ax 


(82) 


V  7c3.45c7?c3.4Tc3.4 

Normal  shock  relations  calculate  the  static  temperature  ratio,  static  pressure  ratio,  total  pres¬ 
sure  ratio,  and  Mach  number  across  the  oblique  shock  to  station  c3.6  [21].  The  entropy  rise  to 
substation  c3.6  is  calculated  using  Eq.  83  [4], 


A5C3.6  =  [Cp,C3.4ln(rc3.6)  —  CP,Rhi(T3.2)\ 


-***(%) 


(83) 


The  speed  of  sound  at  substation  c3.6  is  calculated  using  Eq.  84. 


®c3.6  —  \J 7c3 . 4 fjc > . 4 T f; 


(84) 


The  circumferential  velocity  in  the  detonation  wave  frame-of-reference  VC3.qx,w  is  calculated  as  the 
product  of  the  Mach  number  across  the  oblique  shock  and  ac3.6  from  Eq.  84  [4] . 

The  static  temperature  is  converted  to  stagnation  in  the  detonation  wave  frame-of-reference  to 
get  the  total  temperature  using  Eq.  85. 


—  Tc3,6  + 


V2 

vc3.6x,w 


-  Vlt  +  V i 


3.2  y 


2  Cp, 


c3.4 


(85) 


The  total  pressure  is  calculated  using  isentropic  relations,  as  seen  in  Eq.  86  [21]. 


-ft,c3.6  —  ir3.fi 


/  Tt,c3.6 
\  Tc3  6 


T'cs.r 

7c3.4-! 


(86) 
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The  entropy  rise  after  the  shock  is  calculated  using  Eq.  87  [4]. 


ASt,C3.6  =  [Cp,c3.4  ln(Tt)C3.6)  —  CP,Rln(T3.2)\ 


Rc3A  In 


Pt,c3.6  'N 


Rr  In 


(87) 


Path  d:  Mixed  Flow. 

The  flow  traveling  along  path  d  is  comprised  of  flow  from  substations  a3.4,  63.6,  and  c3.6,  as 
seen  in  Fig.  7.  When  the  three  different  paths  first  meet,  the  parameter  Zdefi  from  Eq.  40  dictates 
what  percentage  of  the  total  flow  is  comprised  of  flow  from  path  c.  Kaemming  assumes  that  the 
remainder  of  the  flow,  calculated  in  Eq.  42,  is  comprised  of  equal  parts  flow  a  and  6  to  simplify  the 
analysis  [4].  Kaemming  uses  the  empirical  models  from  Fig.  9  to  calculate  the  percentage  of  the 
total  flow  coming  from  paths  a  and  6  using  Eq.  88  and  Eq.  89,  respectively  [4].  The  remainder  of 
the  flow  is  comprised  of  flow  from  path  d,  as  seen  in  Eq.  90. 


JJla 

Za  =  — —  =  (1  —  Detonation  Plus  Shock  Flow  Fraction) Zdet  ~  0. 5 Zdefi Relative  Mixing  (88) 

mtot 

ffib 

Zb  =  — - =  Detonation  Plus  Shock  Flow  FractionZ^  —  0.5Zrfej^Relative  Mixing  (89) 

mtot 

Zd  =  ^  =  l-Za-Zb  (90) 

mtot 

Kaemming  uses  the  percentages  from  Eqs.  88-90  to  calculate  weighted  average  values  at  substa¬ 
tion  d3.6  for  the  ratio  of  specific  heat,  specific  gas  constant,  heat  capacity,  total  temperature,  and 
net  entropy  change  using  the  respective  properties  from  stations  a3.4,  63.6,  and  c3.6  [4].  Kaemming 
uses  the  weighted  entropy  change  to  calculate  the  total  pressure  for  the  flow  expanded  to  station 
d3.6  based  on  Eq.  91  [4], 


[Cp,d3. 6  ln(Tt,d3. 6)  — Cp.fl  ln('D.2)+^i?4n(;^p)-AS[j3.6j /Rd3.6| 


(91) 


The  flow  traveling  along  path  d  continues  to  the  exit  of  the  RDE  chamber  with  substation  d3.6 
properties,  as  seen  in  Fig.  7.  After  exiting  the  RDE,  the  flow  is  expanded  to  the  throat  (station  8), 
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exit  (station  9),  and  ambient  using  isentropic  relations  (see  expansion  of  path  a  flow  for  process) 

[4]- 


Exit  Flow  Model 


Kaemming  calculates  the  exit  properties  for  the  three  exiting  paths  a,  b,  and  d  using  weighted 
averages.  The  overall  fuel-to-air  ratio  is  calculated  with  percentages  from  Eq.  40  and  Eq.  42  using 
Eq.  92. 

FARoverau  =  ZdefiFARdefi  +  ZdetFARdet  (92) 

The  specific  thrust  at  the  throat  (^)  8,  the  specific  thrust  at  the  exit  )  9,  the  specific  thrust 

expanded  to  the  ambient  (^)  amb,  static  pressure  at  the  throat  Pg,avg,  static  temperature  at  the 
throat  Ts,avg,  and  the  total  temperature  TtA.avg  are  calculated  using  weighted  averages  with  the 
percentages  from  Eqs.  88-90.  The  specific  impulse  is  calculated  at  the  throat,  exit,  and  ambient 
based  on  Eq.  93  with  i  replaced  with  8,  9,  and  ambient,  respectively  [3;  14]. 


/  F  \  / 1  +  F ARoverau  \  /  gc  \ 

V  rrl  /  avg,i  V  F AR, overall  J  \90  J 


The  equivalent  exhaust  velocity  is  calculated  using  Eq.  94  [3;  14;  22], 


(93) 


Ve  = 


9c 


avg,amb 


(94) 


The  equivalent  exhaust  Mach  number  based  on  the  velocity  from  Eq.  94  is  calculated  using  Eq.  95 


[4]- 

Me  =  /  K  ,  2  =  (95) 

\J lPgORpTtA,avg  ~  (r h) aVg,amb 

The  total  exit  pressure  is  calculated  using  isentropic  relations  with  the  equivalent  exhaust  Mach 
number  from  Eq.  95  [4]. 
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Mass  Flow. 


All  of  the  flow  path  calculations  are  based  on  percentage  of  total  flow.  The  throat  area  is 
based  the  the  cross-sectional  area  from  Eq.  14,  the  user-defined  ratio  of  from  Table  6,  and 

flow  coefficient  for  combustion  chamber  exit  Cws  from  Table  6  [4],  The  exit  mass  flow  (“w4”  in 
Fig.  11)  is  calculated  using  the  throat  area  in  the  function  find_flow  and  correcting  the  result 
with  the  mass  correction  factor  from  Eq.  9.  The  difference  between  the  injection  (“w3”)  and  exit 
mass  (“w4”)  flow  rate  is  updated  and  the  iteration  cycle  varies  the  pre-detonation  pressure  P3.2  to 
equate  the  two  mass  flow  calculations  [4], 

Summary 

The  NPSS  RDE  model  is  composed  of  four  elements:  an  air  source,  a  fuel  source,  the  RDE,  and 
the  exhaust.  The  flow  in  the  RDE  element  itself  is  divided  into  four  paths:  detonation,  detonation 
and  shock,  deflagration,  and  mixed.  The  pre-detonation  pressure  P3.2  is  iterated  until  the  entrance 
and  exit  model  mass  flow  calculations  converge  to  a  consistent  solution. 
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IV.  Analysis  and  Results 


The  NPSS  model  is  verified  by  comparison  to  the  Kaemming  model  [4]  and  published  CFD 
results  [3;  14].  The  default  values  from  Table  6  are  the  baseline  case.  The  baseline  case  is  the 
same  as  the  CFD  simulations  [3;  14],  enabling  a  direct  comparison  between  the  models.  In  all 
of  the  following  cases,  the  default  values  from  Table  6  are  used  for  all  parameters  except  for  the 
parameter  being  varied. 

The  specific  impulses  in  all  of  the  following  cases  is  calculated  for  air-breathing  engines  using 
the  weight  of  the  fuel  alone.  The  specific  impulse  for  rocket  engines  will  have  the  same  magnitude 
as  the  specific  thrust,  but  with  units  of  seconds.  The  use  of  both  specific  impulse  calculations 
allows  the  RDE  to  be  compared  to  both  air-breathing  engines  and  rocket  engines. 

Comparison  to  Kaemming  Model 

The  average  specific  impulse  expanded  to  the  ambient  from  Eq.  93  is  used  as  the  metric  of 
comparison  because  small  changes  are  amplified  due  to  the  1+fARoveraii  term  ( 1  VflR°tiera“  k,  35 
in  the  default  case).  Since  both  the  Kaemming  and  NPSS  models  were  developed  similarly,  the 
results  were  expected  to  be  the  same.  The  only  differences  would  arise  from  any  inherent  numerical 
precision  differences  between  Microsoft  Excel  and  NPSS. 

Figures  12-15  illustrate  four  parametric  variations  to  demonstrate  agreement  between  the  two 
models.  Figures  12  and  13  vary  the  geometry  of  the  RDE  via  the  height  and  the  mean  throughflow 
diameter,  respectively.  Figures  14  and  15  vary  the  thermodynamic  inputs  of  the  RDE  via  the 
entrance  pressure  (Pa)  of  the  air  and  fuel  and  the  entrance  temperature  (Th)  of  the  air  and  fuel, 
respectively.  In  the  analysis,  the  entrance  pressure  and  temperature  for  the  air  and  fuel  are  both 
assumed  to  be  equal  (i.e.,  Pt3,air  =  Pa, fuel  and  Tt3,air  =  Tt3juet )  to  simplify  the  analysis.  As 
expected,  the  results  of  the  Kaemming  model  matched  the  results  of  the  NPSS  model. 
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Figure  12.  Parametric  variation  of  RDE  height  and  specific  impulse  ( Dm  =  140  mm,  Pta 
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Figure  13.  Parametric  variation  of  mean  throughfiow  diameter  and  specific  impulse  (Iirde 
237  mm,  P*3  =  58.784  psia,  and  rt3  =  540  R ) 
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Figure  15.  Parametric  variation  of  entrance  temperature  and  specific  impulse  {Hrde  =  237  mm, 
Drn  =  140  mm,  and  Pt3  =  58.784  psia) 
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Comparison  to  CFD  Simulations 


The  NPSS  model  was  compared  to  published  CFD  results  [3;  14] .  The  parametric  variation  of 
RDE  height  from  Fig.  12  is  compared  to  CFD  results  [14:  Fig.  10]  in  Fig.  16. 


RDE  Height  hRDE  (mm) 

Figure  16.  Comparison  of  the  RDE  height  variation  from  NPSS  model  to  published  results 

( Dm  =  140  mm,  Pt3  =  58.784  psia,  and  Tt3  =  540  R ) 

The  specific  impulse  calculated  from  the  NPSS  model  is  seemingly  invariant  with  RDE  height, 
while  the  CFD  simulation  shows  an  inverse  relationship  between  specific  impulse  and  RDE  height. 
Kaemming  suggests  higher  order  interactions  (e.g.,  non-linearity  of  the  empirical  models  from 
Fig.  9,  calorically  imperfect  gas  effects,  diffusion)  take  place  in  the  CFD  simulations  which  lower 
the  specific  impulse  [4],  Nonetheless,  the  greatest  deviation  between  the  NPSS  model  and  the  CFD 
simulation  is  330  sec,  which  represents  a  9.7%  difference.  While  the  NPSS  model  does  not  follow 
the  CFD  exactly,  it  provides  a  sufficient  approximation. 

The  parametric  variation  of  the  mean  throughflow  diameter  from  Fig.  13  uses  an  entrance  pres¬ 
sure  of  58.784  psia,  as  seen  in  Table  6.  However,  CFD  simulation  results  for  the  mean  throughflow 
diameter  use  an  entrance  pressure  of  147  psia  [14].  The  NPSS  parametric  variation  of  the  mean 
throughflow  diameter  in  Fig.  17  also  uses  an  entrance  pressure  of  147  psia  and  is  compared  to  CFD 
results  [14:  Fig.  8]. 

The  specific  impulses  calculated  from  the  variation  of  mean  diameter  in  the  NPSS  model  agree 
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Figure  17.  Comparison  of  the  RDE  mean  throughflow  diameter  variation  from  NPSS  model 
to  published  results  (Hroe  =  237  mm,  Pt 3  =  147  psia,  and  Tt3  =  540  R) 


fairly  well  with  the  CFD  simulations.  The  greatest  difference  of  280  sec  represents  a  5.5%  difference. 

The  ratio  of  the  entrance  pressure  to  the  ambient  pressure  (Pa/  Pq)  is  used  to  compare  the 
results  from  the  NPSS  model  to  CFD  simulation  [3:  Table  2]  in  Fig.  18.  In  both  the  NPSS  model 
and  CFD  simulations,  the  ambient  pressure  is  Pq  =  14.7  psia  [3;  4],  Besides  the  lowest  pressure 


Pressure  Ratio  P^/Pq 

Figure  18.  Comparison  of  Pa/ Pq  from  NPSS  model  to  published  results  (Iirde  =  237  mm, 
Drn  =  140  mm,  Pq  =  14.7  psia,  and  Tt 3  =  540  R ) 


ratio  of  Ptz/Po  =  2.5, 


the  NPSS  model  consistently  calculates  a  higher  specific  impulse  than  the 
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CFD  simulation.  However,  the  greatest  difference  of  380  sec  represents  a  7.0%  difference,  which 
represents  a  fair  approximation. 

The  NPSS  model  does  not  provide  the  complex  flowfield  visualizations  and  distribution  profiles 
for  properties  of  interest  like  the  CFD  simulations  seen  in  Fig.  2  and  Fig.  3.  However,  the  NPSS 
model  executes  significantly  faster  (less  than  5  sec )  and  requires  fewer  computational  resources 
than  the  CFD  simulations  while  providing  a  sufficient  approximation. 


Model  Application 


The  NPSS  model  was  found  to  agree  with  the  Kaemming  model  [4]  and  published  CFD  sim¬ 
ulation  results  [3;  14].  The  NPSS  model  is  then  used  to  illustrate  why  an  iterative  method  was 
required  to  determine  P3.2  as  well  as  the  effects  of  the  shock  angle  0shock ■  Finally,  the  NPSS  model 
is  used  to  further  analyze  the  RDE  in  a  rocket  configuration. 

As  previously  mentioned,  the  pre-combustion  pressure  P3.2  is  iterated  in  the  RDE  model.  A 
parametric  variation  of  the  entrance  total  pressure  and  entrance  total  temperature  with  the  pre¬ 
combustion  pressure  is  illustrated  in  Fig.  19. 
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Figure  19.  Parametric  variation  of  entrance  total  pressure  and  entrance  total  temperature 
with  the  pre-combustion  pressure  (Jirde  =  237  mm  and  Dm  =  140  mm) 


As  seen  in  in  Fig.  19,  the  relationship  between  the  entrance  total  pressure,  entrance  total 
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temperature,  and  pre-combustion  pressure  is  complex.  At  a  total  entrance  temperature  of  300  R, 
the  relationship  between  the  entrance  total  pressure  and  pre-combustion  pressure  is  non-linear, 
while  at  a  total  entrance  temperature  of  1050  R ,  the  relationship  is  approximately  linear.  Due 
to  the  complex  relationship  between  the  three  parameters,  an  empirical  model  has  not  yet  been 
developed  and  the  iterative  process  is  required. 

Kaemming  uses  a  shock  angle  of  0shock  =  60°  as  the  angle  of  the  shock  relative  to  the  flow  [4], 
as  described  in  Table  6.  A  parametric  analysis  of  the  shock  angle  with  the  total  pressure  and  total 
temperature  exiting  the  RDE  is  illustrated  in  Fig.  20.  A  parametric  analysis  of  the  shock  angle 
with  the  specific  thrust  and  specific  impulse  is  illustrated  in  Fig.  21. 
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Figure  20.  Parametric  analysis  of  shock  angle  with  the  exit  total  pressure  and  total  tempera¬ 
ture  (Jird e  =  237  mm ,  Dm  =  140  mm,  Pt 3  =  58.784  psia,  and  Tt 3  =  540  R) 

As  seen  in  Fig.  20,  the  shock  angle  9shock  does  not  significantly  affect  the  total  pressure  and  total 
temperature  exiting  the  RDE.  The  total  pressure  changes  less  than  1.5%  and  the  total  temperature 
changes  less  than  0.003%  across  a  wide  range  of  shock  angles.  As  seen  in  Fig.  21,  the  shock  angle 
9 shock  does  not  significantly  affect  the  specific  thrust  or  specific  impulse.  The  specific  thrust  changes 
less  than  1%  and  the  specific  impulse  changes  less  than  1.5%  across  a  wide  range  of  shock  angles. 
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Figure  21.  Parametric  analysis  of  shock  angle  with  the  specific  thrust  and  specific  impulse 

(hRDE  =  237  mm,  Dm  =  140  mm,  Pt 3  =  58.784  psia,  and  Tt3  =  540  i?) 


Therefore,  if  there  is  a  deviation  from  from  the  assumed  60°,  the  effect  is  negligible. 

The  effect  of  RDE  geometry  (specifically,  the  RDE  height  and  mean  throughflow  diameter)  on 
the  specific  thrust  and  specific  impulse  is  illustrated  in  Fig.  22  and  Fig.  23,  respectively. 

As  with  the  parametric  analysis  of  the  RDE  height  in  Fig.  12,  the  specific  thrust  in  Fig.  22  and 
the  specific  impulse  in  Fig.  23  reach  plateaus.  Additional  RDE  height  will  not  provide  any  more 
specific  thrust  or  specific  impulse,  so  for  weight  savings,  only  the  shortest  required  RDE  should  be 
used  in  a  rocket.  In  the  plateaus,  the  RDE  has  extracted  as  much  energy  from  the  fuel  as  possible, 
so  additional  time  in  the  combustion  chamber  will  not  yield  additional  energy  extraction.  As  the 
mean  throughflow  diameter  increases,  a  greater  RDE  height  is  required  to  reach  the  plateau.  As 
the  mean  throughflow  diameter  increases,  the  mass  flow  passing  through  the  RDE  also  increases. 
The  greater  mass  flow  requires  additional  time  in  the  combustion  chamber  for  the  RDE  to  extract 
maximum  energy.  Additionally,  as  the  mean  throughflow  diameter  increases,  the  plateau  values 
increase.  It  follows  that  as  the  mass  flow  increases,  the  total  amount  of  energy  extracted  increases. 
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Figure  22.  Parametric  analysis  of  RDE  height  and  mean  throughflow  diameter  on  specific 
thrust  (Pt3  =  58.784  psia ,  and  Tt 3  =  540  R ) 
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Figure  23.  Parametric  analysis  of  RDE  height  and  mean  throughflow  diameter  on  specific 
impulse  (Pts  =  58.784  psia ,  and  Tt 3  =  540  R ) 
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Figure  24  illustrates  the  relationship  between  mean  throughflow  diameter,  RDE  height,  and 
mass  flow  rate.  As  seen  in  Fig.  24,  the  mass  flow  rate  is  invariant  with  the  RDE  height  and  has 


Figure  24.  Parametric  analysis  of  RDE  height  and  mean  throughflow  diameter  on  mass  flow 
rate  (Pt3  =  58.784  psia ,  and  Tt3  =  540  R ) 


a  linear  relationship  with  the  mean  throughflow  diameter.  The  linear  relationship  between  mean 
throughflow  diameter  and  mass  flow  rate  provides  a  simple  scaling  model  that  can  be  used  to 
determine  the  desired  size  of  an  RDE. 

Figure  25  illustrates  the  relationship  between  mean  throughflow  diameter,  RDE  height,  and 
thrust.  After  a  certain  height,  the  specific  thrust  in  Fig.  22  and  the  thrust  in  Fig.  25  both  reach 


Figure  25.  Parametric  analysis  of  RDE  height  and  mean  throughflow  diameter  on  thrust 

(Pt 3  =  58.784  psia,  and  Tt3  =  540  R) 


plateaus.  In  designing  an  RDE,  only  the  minimum  required  RDE  height  (i.e.,  the  height  where  the 


48 


specific  thrust  and  thrust  reach  plateaus)  should  be  used;  there  are  no  additional  gains  from  a  taller 
RDE.  In  Fig.  22,  increasing  the  mean  throughflow  diameter  results  in  decreasing  gains  in  specific 
thrust.  In  Fig.  25,  however,  increasing  the  mean  throughflow  diameter  results  in  increasing  gains 
in  the  thrust.  Based  on  Figs.  22-25,  the  mean  throughflow  diameter  is  the  key  sizing  parameter. 

The  combined  effect  of  total  entrance  temperature  and  total  entrance  pressure  on  the  specific 
thrust  and  specific  impulse  is  illustrated  in  Fig.  26  and  Fig.  27,  respectively. 


Figure  26.  Parametric  analysis  of  total  entrance  temperature  and  total  pressure  on  specific 
thrust  (Iirde  =  237  mm  and  Dm  =  140  mm) 

At  the  lower  entrance  pressures,  the  total  entrance  temperature  has  a  significant  effect  on  the 
specific  thrust  and  specific  impulse.  As  the  entrance  total  pressure  increases,  the  total  entrance 
temperature  has  less  of  an  effect.  Increasing  the  total  entrance  temperature  decreases  the  specific 
thrust  and  specific  impulse  for  the  different  total  entrance  pressures.  The  specific  thrust  and  specific 
impulse  are  increasing  at  a  decreasing  rate  and  appear  to  be  trending  towards  plateaus  at  150  ,, 

I'Om/  S 

and  5500  sec,  respectively. 
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Figure  27.  Parametric  analysis  of  total  entrance  temperature  and  total  pressure  on  specific 
impulse  (Iirde  =  237  mm  and  Dm  =  140  mm) 


The  combined  effect  total  entrance  temperature  and  total  entrance  pressure  on  the  mass  flow  is 
illustrated  in  Fig.  28.  The  mass  flow  rate  from  Fig.  28  is  used  with  the  specific  thrust  from  Fig.  26 


Figure  28.  Parametric  analysis  of  total  entrance  temperature  and  total  pressure  on  mass  flow 
rate  (Iirde  =  237  mm  and  Dm  =  140  mm) 


to  calculate  the  total  thrust,  as  seen  in  Fig.  29. 

Both  the  mass  flow  rate  in  Fig.  28  and  thrust  in  Fig.  29  increase  linearly  with  the  total  entrance 
pressure.  The  entrance  total  temperature  affects  the  slope  of  each  line,  and  at  a  higher  total  entrance 
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Figure  29.  Parametric  analysis  of  total  entrance  temperature  and  total  pressure  on  thrust 

(h.RDE  =  237  mm  and  Dm  =  140  mm) 

pressure,  the  total  entrance  temperature  has  a  greater  impact.  Even  if  the  specific  thrust  plateaus 
as  in  Fig.  26,  the  increasing  mass  flow  rate  will  increase  the  net  thrust. 

Material  properties  often  limit  the  temperature  exiting  the  combustor  in  traditional  engines 
[7].  Therefore,  the  temperature  exiting  the  RDE  7*4  is  a  parameter  of  concern.  The  effect  of  total 
entrance  temperature  and  total  entrance  pressure  Pt3  on  total  exit  temperature  Tf4  is  illustrated 
in  Fig.  30. 


Figure  30.  Parametric  analysis  of  total  entrance  temperature  and  total  pressure  on  total  exit 
temperature  {Jird e  =  237  mm  and  Dm  =  140  mm) 
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The  total  entrance  pressure  has  a  negligible  effect  on  the  total  exit  temperature.  There  is  a  linear 
relationship  between  the  total  entrance  temperature  and  total  exit  temperature.  Temperatures  of 
3000  R  to  4700  R  exist  in  real  rocket  engines  [22],  though  material  developments  may  lead  to  higher 
allowable  temperatures.  The  total  exit  temperatures  in  Fig.  30  fall  within  the  range  of  existing 
rocket  exhaust  temperatures;  however,  cooling  will  still  be  a  concern. 

The  RDE  is  a  pressure-gaining  combustor,  which  is  why  the  RDE  is  able  to  achieve  the  higher 
efficiency  of  the  Atkinson  cycle  [3;  11;  13;  14].  The  pressure  ratio  across  the  RDE  is  given  in 
Eq.  96. 


Pi 


-Kb  = 


t.  8 


t,  3.2 


(96) 


The  parametric  variation  of  Kb  with  total  entrance  pressure  and  total  entrance  temperature  is 
illustrated  in  Fig.  31.  The  value  of  Kb  =  1  is  shown  to  illustrate  that  all  pressure  ratios  are  greater 
than  one. 


Figure  31.  Parametric  variation  of  Kb  with  total  entrance  pressure  and  total  entrance  temper¬ 
ature  (Iirde  =  237  mm  and  Dm  =  140  mm) 


As  seen  in  Fig.  31,  the  RDE  produces  a  pressure  rise  across  the  entire  range  of  pressures  and 
temperatures  implemented.  The  pressure  rise  agrees  with  CFD  simulations  [11;  14].  Therefore, 
the  RDE  does,  in  fact,  follow  the  Atkinson  cycle  and  will  be  able  to  realize  the  higher  thermal 
efficiency  [3;  11;  13;  14]. 

The  pressure  profile  from  Fig.  8  results  in  a  pressure  loss  across  the  plenum  from  station  3  to 
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Figure  32.  Parametric  variation  of  plenum  pressure  loss  with  total  entrance  pressure  and  total 
entrance  temperature  {Hrde  =  237  mm  and  Dm  =  140  mm) 

substation  3.2.  The  parametric  variation  of  the  plenum  pressure  loss  with  total  entrance  pressure 
and  total  entrance  temperature  is  illustrated  in  Fig.  32. 

As  seen  in  Fig.  32,  the  pressure  loss  is  significant;  the  total  pressure  that  enters  the  RDE 
combustion  chamber  (substation  3.2)  is  only  25%  to  35%  of  the  total  pressure  entering  the  RDE 
(station  3).  Some  of  the  pressure  loss  can  be  attributed  to  feedback  from  the  injectors  [17].  Detailed 
analysis  of  the  pressure  loss  across  the  plenum  requires  further  investigation  and  is  outside  the  scope 
of  this  thesis. 


Investigative  Questions  Answered 

Ibf 

The  highest  specific  thrust  calculated  is  just  over  150  T, — f—  and  the  highest  specific  impulse 
is  approximately  5,500  sec ,  as  seen  in  Fig.  26  and  Fig.  27,  respectively.  Again,  the  specific  impulse 
is  force  per  weight  flow  of  fuel  alone.  The  specific  impulse  based  on  total  air  and  fuel  weight  flow, 
for  comparison  to  other  rocket  propulsion  methods,  is  150  sec. 

The  specific  impulse  for  a  traditional  turbojet  using  hydrogen  fuel  can  range  from  4,500  sec 
to  7,000  sec  [23].  The  calculated  RDE  specific  impulse  for  air-breathing  engines  of  5,500  sec  falls 
within  the  range  of  turbojet  engines.  A  typical  rocket  engine  using  liquid  hydrogen  fuel  has  a 
specific  impulse  of  410  sec  [22:  Table  1.6].  The  calculated  RDE  specific  impulse  for  rocket  engines 
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of  150  sec  falls  below  that  of  existing  rocket  engines.  The  nozzle  and  exhaust  were  assumed  to  be 
ideal,  as  were  the  internal  flows  between  the  air  and  fuel  sources  into  the  RDE.  The  RDE  requires 
further  investigation  in  order  to  improve  the  performance  of  the  RDE.  In  a  real  RDE,  the  specific 
impulse  will  be  lower  due  to  losses,  but  the  results  are  still  very  promising  and  merit  continued 
development  of  the  RDE  in  a  rocket  configuration. 

Summary 

The  NPSS  model  was  verified  by  comparison  to  the  Kaemming  model  [4]  and  published  CFD 
results  [3;  14].  The  NPSS  model  matched  the  Kaemming  model.  The  NPSS  model  matched  CFD 
simulations  closely  with  a  maximum  difference  of  10%.  The  calculated  specific  impulses  of  the  RDE 
peaked  at  5,500  sec  when  based  on  the  weight  flow  of  fuel  alone  and  150  -sec  when  based  on  the 
weight  flow  of  fuel  and  air.  The  promising  results  merit  continued  development  of  the  RDE  model 
and  development  of  a  real  RDE  in  a  rocket  configuration  for  comparison. 
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V.  Conclusions  and  Recommendations 


Detonation-based  combustion  follows  an  Atkinson  cycle  and  has  a  higher  inherent  thermody¬ 
namic  efficiency  than  deflagration-based  combustion  following  a  Brayton  cycle  [8].  The  parametric 
variations  in  this  thesis  calculated  that  the  RDE  had  an  air-breathing  engine  specific  impulse  that 
fell  in  the  range  of  traditional  turbojets  using  hydrogen  fuel  and  a  rocket  engine  specific  impulse 
that  fell  below  hydrogen-oxygen  rockets.  Nonetheless,  the  advantages  of  the  RDE  are  numerous. 
For  one,  a  higher  specific  impulse  would  require  a  lower  mass  flow  to  produce  the  same  amount  of 
thrust  [22],  and  lower  mass  flow  translates  to  a  lower  total  required  fuel  load. 

Recommendations  for  Future  Research 

There  are  two  main  categories  for  recommendations:  model  development  and  model  applica¬ 
tion.  Model  development  describes  improvements  that  can  be  made  to  the  NPSS  model.  Model 
application  describes  possible  uses  for  the  NPSS  model. 

Model  Development. 

While  the  simplified  NPSS  model  agrees  with  the  CFD  simulations,  it  has  not  been  compared 
to  experimental  results.  The  empirical  models  are  based  exclusively  on  results  from  higher-order 
CFD  simulations  [3;  14].  Additionally,  the  exit  flow  boundary  conditions  in  the  CFD  simulations 
are  still  a  matter  of  debate  [3;  4;  14].  The  empirical  models,  and  thus  the  entire  RDE  model,  would 
be  further  improved  if  experimental  results  were  incorporated  into  their  development.  Including 
experimental  results  would  also  validate  the  code. 

The  NPSS  model  is  programmed  using  iterative  loops.  The  mass  flows  for  air  and  fuel  are 
calculated  in  those  iterative  loops,  as  are  the  fuel-to-air  ratios  for  the  deflagrated  and  detonated 
flows.  Additionally,  the  inlet  flow  conditions  ( e.g pressure,  temperature)  are  coded  into  RDE 
element.  The  NPSS  environment  allows  for  all  those  parameters  to  be  defined  outside  the  RDE 
element  in  the  air  and  fuel  inlet  elements.  Also,  NPSS  contains  solvers  to  facilitate  the  iteration 
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loops.  However,  the  current  code  does  not  incorporate  these  capabilities.  Future  improvements  to 
the  code  should  incorporate  the  NPSS  capabilities  to  simplify  the  code  and  speed  up  calculations. 

Finally,  the  RDE  model  assumes  the  flow  is  calorically  perfect  and  has  constant  chemical  com¬ 
position  [4],  In  traditional  deflagration  combustion  cycle  analysis,  equilibrium  chemistry  is  often 
assumed  for  a  more  accurate  calculation  [4;  18] .  Equilibrium  chemistry  should  be  an  incorporated 
capability  of  the  RDE  model. 

Model  Application. 

The  rocket  configuration  outlined  in  this  thesis  is  a  starting  point  for  the  implementation  of 
the  RDE.  Further  analysis  should  be  accomplished  on  the  different  configuration  options.  Several 
configurations  utilizing  a  generic  combustor  have  been  developed  [16].  The  RDE  model  can  be 
used  in  conjunction  with  a  compressor  to  increase  the  pressure  going  into  the  RDE.  Another 
option  is  to  couple  the  RDE  with  a  traditional  rocket  engine  so  the  RDE  can  be  used  in  parallel 
with  the  traditional  combustor  to  create  additional  thrust.  Alternatively,  it  can  be  used  in  an 
afterburner- like  capacity  to  further  increase  thrust. 

A  final  potential  use  for  the  RDE  is  in  a  traditional  turbojet  engine.  It  can  be  used  to  replace 
the  burner  and  compressor  in  the  engine,  or  just  the  burner  itself.  Alternatively,  it  could  be  used 
as  a  secondary  combustor,  either  in  parallel  with  a  traditional  combustor  in  an  afterburner-like 
capacity. 

The  NPSS  RDE  element  allows  the  various  configurations  to  be  readily  modeled.  The  analysis 
could  be  used  to  find  optimal  design  and  operating  conditions  for  the  RDE.  Additionally,  the 
analyses  would  provide  a  more  direct  comparison  between  traditional  engines  and  those  modified 
with  an  RDE. 

Summary 

The  RDE  provides  promising  results  that  merit  further  development.  The  implemented  empir¬ 
ical  models  and  gas  dynamic  properties  can  be  improved  with  experimental  results  to  provide  more 
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accurate  results.  The  internal  structure  of  the  RDE  code  can  be  improved  and  simplified  by  using 
some  of  the  built-in  NPSS  features.  The  NPSS  RDE  element  allows  the  various  configurations  to 
be  readily  modeled  to  find  optimal  design  and  operating  conditions  for  their  various  uses. 
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Appendix  A.  Summary  of  Implemented  Functions 


Several  functions  are  implemented  for  repeated  calculations.  Functions  are  coded  in  the  same 
manner  to  subroutines.  The  key  difference  is  that  functions  return  values  and  subroutines  do  not; 
in  the  RDE  element,  all  functions  are  of  type  real,  while  all  subroutines  are  of  type  void.  The 
functions  were  used  to  perform  common  thermodynamic  calculations. 


Isentropic  Relations 


The  functions  PtqP  and  TtqT  solved  for  the  isentropic  total/static  pressure  and  temperature  ra¬ 
tios,  respectively.  The  isentropic  total/static  pressure  was  calculated  based  on  Eq.  A.l  [21:  Eq.  (44)]. 


A 

P 


(A.l) 


The  isentropic  total/static  temperature  was  calculated  based  on  Eq.  A. 2  [21:  Eq.  (43)]. 


Tt 

T 


=  1  + 


(A.2) 


The  functions  M_PtqP  and  M_TtqT  solved  for  the  Mach  number  based  on  the  isentropic  to¬ 
tal/static  pressure  and  temperature  ratios,  respectively,  based  on  Eqs.  A. 3  and  A. 4,  respectively. 


M  = 


(A.3) 


M  = 


(A.4) 


The 
Eq.  A. 5 


function  AqAstar  solved  for  the  isentropic  area  ratio  for  chocked  flow  in  a  nozzle,  based  on 
[21:  Eq.  (80)]. 


A 

A* 


(A.5) 
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The  functions  Msub_A  and  Msup_A  utilize  the  function  AqAstar  to  iteratively  calculate  the  Mach 
number  based  on  the  isentropic  area  ratio  and  produce  the  subsonic  and  supersonic  solution,  re¬ 
spectively. 


Normal  Shock  Relations 


The  functions  P2qPlns,  Pt2qPtlns,  T2qTlns,  and  M2ns  calculated  flow  properties  across  a 
normal  shock.  The  static  pressure  ratio  across  a  normal  shock  was  calculated  by  the  function 
P2qPlns  using  Eq.  A. 6  [21:  Eq.  (93)]. 

ft  =  27M2  -  (7  -  1)  (A. 6) 

Pi  7  +  1 


The  total  pressure  ratio  across  a  normal  shock  was  calculated  by  the  function  Pt2qPtlns  using 
Eq.  A. 7  [21:  Eq.  (99)]. 


Pt2 

(7  +  1)M2 

"Y 

7-1 

7+l 

Pn  ~ 

_(7-  l)M2  +  2_ 

2  7M2  —  (7  —  1) 

(A.7) 


The  static  temperature  ratio  across  a  normal  shock  was  calculated  by  the  function  T2qTlns  using 
Eq.  A. 8  [21:  Eq.  (95)]. 

T2  _  [27M2  -  (7  -  1)][(7  -  m2  +  2]  /  *  q\ 

Ti  (7+l)2M2  1  j 


The  Mach  number  across  a  normal  shock  was  calculated  by  the  function  M2ns  using  Eq.  A. 9 


[21:  Eq.  (96)]. 


Ah 


/ 


(7  -  1  )Ml  +  2 

27M2-(7-l) 


(A.9) 


Mass  Flow 


The  function  f  ind_f  low  iterates  total  mass  flow  to  match  a  desired  total  area  for  three  separate 
flows.  The  function  uses  the  desired  flow  area  and  three  arrays  containing  the  density,  velocity,  and 
flow  percentages  of  each  of  the  three  streams.  The  area  contribution  each  of  the  flows  made  was 
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calculated  using  Eq.  A.  10. 


(A.10) 


Ai 


mtot(%  flow)j 
pAA% 


The  total  mass  flow  was  iterated  until  the  desired  total  area  At0t  was  attained  by  Eq.  A.ll. 


Atot  =  Ai  +  A2  +  A3 


(A.ll) 


Nozzle  Expansion 


The  function  expans ionunode  finds  the  exit  Mach  number  for  a  flow  expanding  from  the  throat 
to  exit  (station  8  to  station  9)  by  first  calculating  the  isentropic  area  ratio  ,  and  then  using 


Eq.  A. 12. 


(A.12) 


The  subsonic  Mach  number  at  station  9  was  calculated  using  Msub_A.  The  subsonic  Mach  number 
was  used  to  calculate  Pg  using  Eq.  A.  13. 


Py  =  Ps 


(A.13) 


The  calculated  pressure  Pg  was  compared  to  the  ambient  pressure:  if  it  was  less  than  or  equal  to 
the  ambient  pressure  Po,  the  expansion  was  subsonic. 

If  the  calculated  pressure  P9  was  greater  than  the  ambient  pressure  Po,  the  expansion  was 
supersonic.  The  supersonic  Mach  number  at  station  9  was  calculated  using  Msup_A,  and  the  pressure 
at  station  9  was  calculated  using  Eq.  A.13.  The  normal  shock  wave  was  accounted  for  by  multiplying 
the  calculated  pressure  by  the  appropriate  factor  for  a  normal  shock  using  the  function  P2qPlns. 
If  the  calculated  pressure  Pg  after  the  normal  shock  was  still  greater  than  the  ambient  pressure  Po, 
the  nozzle  was  fully  supersonic. 

If  the  calculated  pressure  Pg  after  the  normal  shock  was  less  than  the  ambient  pressure  Po, 
the  nozzle  was  overexpanded  and  there  was  a  shock  in  the  nozzle.  The  equation  used  to  find  the 
exit  Mach  number  was  derived  from  the  continuity  equation  (i.e.,  ms  =  mg),  based  on  Eq.  A. 14 
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[7:  Eq.  (1.4)]. 


m  =  PAM 


(A.14) 


Assuming  the  ratio  of  specific  heats  and  the  gas  constant  remain  constant  during  the  expansion, 
the  full  continuity  equation  becomes  Eq.  A.  15. 


PgAgMs 


PgAgMg 


(A.15) 


The  pressure  at  the  exit  was  equal  to  the  ambient  pressure  ( Pg  =  Po).  Assuming  chocked  flow  at 
the  nozzle  (Ms  =  1)  and  adiabatic  expansion  (Tts  =  To),  the  equation  reduces  to  Eq.  A.  16. 


TsA8 


/ 


PgAgMg 


(A.16) 


Squaring  both  sides  and  rearranging  leads  to  Eq.  A. 17,  which  was  used  to  iteratively  solve  for  Mg. 


Mi 


(A.17) 


Chapman-Jouguet  with  Area  Relief 

The  function  CJjiach_area  found  the  Chapman-Jouguet  Mach  number  for  a  flow  with  two 
different  values  of  7  and  a  changing  area.  A  steady,  one-dimensional  detonation  propagates  at  a 
sonic  velocity  at  the  Chapman-Jouguet  state  [19].  In  the  rotating  detonation  engine,  however,  there 
is  lateral  relief  after  the  detonation  which  provides  another  source  of  energy  relief  [4].  The  lateral 
relief  mode  was  modeled  by  an  increase  in  area  after  the  detonation,  which  produced  the  slower 
detonation  and  decreased  pressure  rise  provided  by  the  lateral  relief  [4]. 

The  mass,  momentum,  and  energy  conservation  equations  were  implemented  with  a  change  in 
area  and  gas  properties  from  state  1  to  2.  The  equations  of  conservation  of  mass,  momentum,  and 
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energy  became  Eqs.  A.  18,  A.  19,  and  A. 20,  respectively  [4], 


Cp\T\ 


P\  A  i  M ! 


P2A2M2 


PiA\(l  +  71  Mi)  =  P2A2 


1  + 


CP2T2 


/  7250 

V  F2T2 

(A.18) 

+  72 -M2  ^ 

(A.19) 

[l+(72  2  )"2] 

(A. 20) 

In  the  CFD  simulations,  the  heat  of  combustion  is  corrected  to  stoichiometric  heat  addition  per 
total  reactant  flow.  The  traditional  heat  of  combustion  AH  of  the  fuel  is  related  to  the  heat  of 
combustion  referenced  to  fuel  flow  A Hv  and  stoichiometric  heat  of  combustion  A Hst0ich  in  Eq.  A. 21 
[4;  15;  18], 


m  fuel  AH v 


(mair  +  rrifuei )  AH  =  (mair  +  m.fuel) 


f  A  H^ioich  \ 

— -  mf 

V  stench  J 


(A.21) 


Dividing  Eq.  A.21  by  the  mass  of  air  results  in  Eq.  A. 22. 


F  AR(ietAHv 


(1  +  F ARdet)  AH  —  (1  +  F ARdet)  AHst(yich 


{  FARdet  \ 

\  F ARgtcieh  J 


(A. 22) 


Kaemming  defines  the  non-dimensional  heat  addition  term  q  in  Eq.  A. 20  as  Eq.  A. 23  using  the 
same  values  for  Cp  1  and  T\  as  in  Eq.  A. 20  [4;  15]. 


_  F ARdetAHy 
CpjT] 


Using  the  definitions  from  Eq.  A. 22,  Eq.  A. 23  can  be  rewritten  using  known  properties,  as  seen  in 
Eq.  A. 24. 


q  = 


AH  (1  + FARdet) 


(A. 24) 


CpjT, 

The  Mach  number  at  station  2  was  fixed  to  one  because  the  detonation  moves  sonically  [4;  18; 
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19].  Rearranging  Eq.  A. 20  to  solve  for  the  temperature  ratio  leads  to  Eq.  A. 25 


^-1 

(  Cpi\ 

"l+  1 

(v) 

1  Ml  +  q 

Ti  1 

\Cp2  J 

72+1 

2 

(A. 25) 


Dividing  Eq.  A.  18  by  Eq.  A.  19  results  in  Eq.  A. 26. 


1  +  7i  M'l 

Wl 


(A. 26) 


Substituting  Eq.  A. 25  into  Eq.  A. 26  and  using  the  relation  ^  =  +L_  [21:  Eq.  (16)]  results  in 
Eq.  A. 27  [4], 


l  +  7i  Mi2 

Mi 


f  72  -  1 
V7i  -  1 


1  +  (^)  M,2  +  q 

72  +  1 
2 


The  function  CJ_mach_area  iterates  Mi  to  satisfy  the  Eq.  A. 27. 


(A. 27) 


Interpolation 

NPSS  allows  for  lookup  and  interpolation  in  tables,  but  does  not  allow  for  arrays  to  be  converted 
into  tables.  The  function  OneDNtrp  was  used  to  interpolate  within  an  array.  The  function  went 
through  each  cell  of  the  independent  array  to  determine  which  two  cells  the  desired  value  was 
between.  Then,  it  interpolated  for  the  desired  value  in  the  corresponding  dependent  array.  The 
function  also  allowed  for  extrapolation;  however,  since  all  of  the  implementations  of  the  function 
were  for  a  bounded  range  for  a  ratio  between  zero  and  one,  no  extrapolation  was  utilized. 


63 


References 


1.  Daniel  Paxson.  RDE  Thermodynamic  Cycle  Model  Status.  Private  communication,  November 

2013. 

2.  Craig  A  Nordeen,  Douglas  Schwer,  Fredrick  Schauer,  John  Hoke,  Baki  Cetegen,  and  Thomas 
Barber.  Thermodynamic  Modeling  of  a  Rotating  Detonation  Engine.  In  49th  AIAA  Aerospace 
Sciences  Meeting,  Orlando,  FL ,  Orlando,  FL,  January  2011.  American  Institute  of  Aeronautics 
and  Astronautics. 

3.  Douglas  A  Schwer  and  Kailas  Kailasanath.  Numerical  Investigation  of  Rotating  Detonation 
Engines.  AIAA  Paper  2010-6880,  2010. 

4.  Tom  Kaemming.  A  User’s  Guide  to  a  Rotating  Detonation  Engine  Thermodynamic  Cycle 
Performance  Code.  Version  Q.  Draft.  Also  refers  to  associated  Microsoft  Excel  code.,  2014. 

5.  Craig  A  Nordeen,  Douglas  Schwer,  Fredrick  Schauer,  John  Hoke,  B  Cetegen,  and  T  Barber. 
Divergence  and  Mixing  in  a  Rotating  Detonation  Engine.  In  51st  AIAA  Aerospace  Sciences 
Meeting,  Grapevie,  TX,  January  2013.  American  Institute  of  Aeronautics  and  Astronautics. 

6.  Z.  Spakovszky.  Thermodynamics  and  Propulsion.  Online,  2007.  Available  at 
http://web.mit.edu/16.unified/www/SPRING/propulsion/notes/node27.html. 

7.  Jack  Mattingly.  Aircraft  Engine  Design.  American  Institute  of  Aeronautics  and  Astronautics, 
Reston,  VA,  2002.  ISBN  9781563475382. 

8.  Daniel  E  Paxson.  A  Simplified  Model  for  Detonation  Based  Pressure-Gain  Combustors.  Na¬ 
tional  Aeronautics  and  Space  Administration,  Glenn  Research  Center,  2010. 

9.  BV  Voitsekhovskii.  Maintained  Detonations.  Soviet  Physics  Doklady,  4,  1960. 

10.  Frank  K  Lu  and  Eric  M  Braun.  Rotating  Detonation  Wave  Propulsion:  Experimental  Chal¬ 
lenges,  Modeling,  and  Engine  Concepts.  Journal  of  Propulsion  and  Power,  30(5):1125-1142, 

2014. 

11.  Craig  A  Nordeen.  Thermodynamics  of  a  Rotating  Detonation  Engine.  PhD  thesis,  University 
of  Connecticut,  2013. 

12.  K  Kailasanath.  The  Rotating-Detonation-Wave  Engine  Concept:  A  Brief  Status  Report.  In 
49th  AIAA  Aerospace  Sci.  Meeting,  Orlando,  FL,  January  2011.  American  Institute  of  Aero¬ 
nautics  and  Astronautics. 

13.  Craig  A  Nordeen,  Douglas  Schwer,  Fredrick  Schauer,  John  Hoke,  Baki  Cetegen,  and  Thomas 
Barber.  Energy  Transfer  in  a  Rotating  Detonation  Engine.  In  J7f/t  AIAA  Joint  Propulsion 
Conference,  San  Diego,  CA,  August  2011.  American  Institute  of  Aeronautics  and  Astronautics. 

14.  Douglas  A  Schwer  and  Kailas  Kailasanath.  Numerical  Study  of  the  Effects  of  Engine  Size 
on  Rotating  Detonation  Engines.  In  49th  AIAA  Aerospace  Sciences  Meeting ,  Orlando,  FL, 
January  2011.  American  Institute  of  Aeronautics  and  Astronautics. 


64 


15.  Tom  Kaemming.  RE:  Cold  Flow.  Private  communication,  January  2015. 

16.  Thomas  A  Kaemming,  Gary  L  Lidstone,  and  Brad  C  Sammann.  Proposed  Nomenclature  Guide 
for  Pulse  Detonation  Engines.  In  38th  AIAA/ASME/SAE/ASEE  Joint  Propulsion  Conference 
and  Exhibit,  Indianapolis,  IN,  July  2002.  American  Institute  of  Aeronautics  and  Astronautics. 

17.  Douglas  A  Schwer  and  Kailas  Kailasanath.  Feedback  into  Mixture  Plenums  in  Rotating  Deto¬ 
nation  Engines.  In  50th  AIAA  Aerospace  Sciences  Meeting  including  the  New  Horizons  Forum 
and  Aerospace  Exposition,  Nashville,  TN,  January  2012.  American  Institute  of  Aeronautics  and 
Astronautics. 

18.  Stephen  Turns.  An  Introduction  to  Combustion:  Concepts  and  Applications.  McGraw-Hill, 
New  York,  3  edition,  2012.  ISBN  9780073380193. 

19.  Paul  W  Cooper.  Explosives  Engineering.  Wiley-VCH,  New  York,  1996. 

20.  Jozef  Jarosinski.  Combustion  Phenomena  Selected  Mechanisms  of  Flame  Formation,  Propaga¬ 
tion,  and  Extinction.  CRC  Press,  Boca  Raton,  2009.  ISBN  9780849384097. 

21.  Ames  Research  Staff.  Equations,  Tables,  and  Charts  for  Compressible  Flow.  Technical  Report 
1135,  National  Advisory  Committee  for  Aeronautics,  1953. 

22.  Jack  Mattingly.  Elements  of  Propulsion:  Gas  Turbines  and  Rockets.  American  Institute  of 
Aeronautics  and  Astronautics,  Reston,  VA,  2006.  ISBN  9781563477799. 

23.  The  University  of  Texas  at  Arlington  Aerodynamic  Research  Center.  Pulsed  Detonation  En¬ 
gines.  Online,  2014.  Available  at 

http://arc.uta.edu/research/pde.htm. 


65 


Vita 


First  Lieutenant  Nihar  N.  Shah  graduated  in  2011  from  Rensselaer  Polytechnic  Institute  in 
Troy,  NY  with  a  Bachelor  of  Science  degree  in  aeronautical  and  mechanical  engineering.  He  was 
commissioned  as  a  second  lieutenant  in  the  United  States  Air  Force  at  Detachment  550  AFROTC 
at  Rensselaer  Polytechnic  Institute.  His  first  assignment  was  to  the  Air  Force  Test  and  Evaluation 
Center  Detachment  5  at  Edwards  AFB  as  an  operational  test  engineer  focusing  on  communication 
system  upgrades  on  bomber  aircraft.  In  August  of  2013,  Lieutenant  Shah  entered  the  Air  Force 
Institute  of  Technology  in  pursuit  of  a  Master  of  Science  degree  in  aeronautical  engineering.  Upon 
graduation,  he  will  be  assigned  to  the  National  Air  and  Space  Intelligence  Center. 


66 


REPORT  DOCUMENTATION  PAGE 


Form  Approved 
OMB  No.  0704-0188 


The  public  reporting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information,  including 
suggestions  for  reducing  this  burden  to  Department  of  Defense,  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports  (0704—0188),  1215  Jefferson  Davis  Highway, 
Suite  1204,  Arlington,  VA  22202—4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  any  penalty  for  failing  to  comply  with  a  collection 
of  information  if  it  does  not  display  a  currently  valid  OMB  control  number.  PLEASE  DO  NOT  RETURN  YOUR  FORM  TO  THE  ABOVE  ADDRESS. 

1.  REPORT  DATE  (DD-MM-YYYY)  2.  REPORT  TYPE  3.  DATES  COVERED  (From  —  To) 

26-03-2015  Master’s  Thesis  Sept  2013  —  Mar  2015 

4.  TITLE  AND  SUBTITLE  5a.  CONTRACT  NUMBER 


Computer  Modeling  of  a  Rotating  Detonation  Engine  in  a  Rocket 

Configuration 


5b.  GRANT  NUMBER 


I  5c.  PROGRAM  ELEMENT  NUMBER 


6.  AUTHOR(S) 


I  5d.  PROJECT  NUMBER 


Shah,  Nihar 


I  5e.  TASK  NUMBER 


I  5f.  WORK  UNIT  NUMBER 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

Air  Force  Institute  of  Technology 

Graduate  School  of  Engineering  and  Management  (AFIT /EN) 
2950  Hobson  Way 
WPAFB  OH  45433-7765 

9.  SPONSORING  /  MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

Air  Force  Research  Lab  Aerospace  Directorate 
Combustion  Branch,  Turbine  Engine  Division 
1790  Loop  Road  North 
WPAFB  OH  45433-7765 


8.  PERFORMING  ORGANIZATION  REPORT 
NUMBER 


AFIT-ENY-MS-15-M-230 


10.  SPONSOR/MONITOR’S  ACRONYM(S) 

AFRL/RQTC 

11.  SPONSOR/MONITOR’S  REPORT 
NUMBER(S) 


12.  DISTRIBUTION  /  AVAILABILITY  STATEMENT 

DISTRIBUTION  STATEMENT  A: 

APPROVED  FOR  PUBLIC  RELEASE;  DISTRIBUTION  UNLIMITED. 

13.  SUPPLEMENTARY  NOTES 

This  work  is  declared  a  work  of  the  U.S.  Government  and  is  not  subject  to  copyright  protection  in  the  United  States. 

14.  ABSTRACT 

Detonation-based  combustors  leverage  the  higher  thermodynamic  efficiency  of  the  Atkinson  cycle  compared  to  the 
traditional  deflagration-based  combustion  of  the  Brayton  cycle.  The  rotating  detonation  engine  (RDE)  has  one  or  more 
shock  waves  rotating  around  an  annulus.  The  RDE  can  theoretically  be  20%  more  thermally  efficient  than  a  traditional 
deflagration-based  cycle.  A  RDE  was  modeled  in  Numerical  Propulsion  System  Simulation  (NPSS)  based  on  a  model 
developed  in  Microsoft  Excel.  The  thermodynamic  analysis  of  the  RDE  in  these  models  is  broken  into  four  streams. 
Empirical  models  were  used  to  find  the  percentage  of  the  total  flow  in  each  stream.  The  pre-detonation  pressure  was 
iterated  until  the  entrance  mass  flow  calculations  matched  the  exit  mass  flow  calculations.  A  parametric  analysis  was 
used  to  compare  the  variation  in  specific  impulse  from  the  NPSS  model  to  the  Microsoft  Excel  model  and  other  published 
results.  The  RDE  has  a  peak  air-breathing  engine  specific  impulse  of  approximately  5,500  sec  and  a  peak  rocket  engine 
specific  impulse  of  approximately  150  sec. 

15.  SUBJECT  TERMS 

Rotating  Detonation  Engine,  Pressure-Gain  Combustor,  Atkinson  Cycle,  Detonation  Waves,  Combustion, 

Thermodynamic  Cycles,  Rocket  Engines 


16.  SECURITY  CLASSIFICATION  OF: 


a.  REPORT 


b.  ABSTRACT!  c.  THIS  PAGE 


17.  LIMITATION  OF 
ABSTRACT 


18.  NUMBER  19a.  NAME  OF  RESPONSIBLE  PERSON 
PAGES  Dr'  R  L  Kin&  AFIT/ENY 

19b.  TELEPHONE  NUMBER  (include  area  code) 

ot)  (937)  785-3636,  x4628;  paul.king@afit.edu 


Standard  Form  298  (Rev.  8-98) 

Prescribed  by  ANSI  Std.  Z39.18 


